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CONVERSION  FACTORS 


lmb  -  1000  dynes/an2  -  .029532  in.  Hg  -  0.1  kilopascal  (kpa) 

2 

1  inch  mercury  (Hg)  •  33.86389mb  _  33863.89  dynes/ an 

2-5  2 

1  dyne  *  1  g  cm/s  -  10  kg  m/s 

1  erg  «  1  dyne  an  =  2.3892  X  10  ®  calorie 

1  foot  =  0.3048m  =  30.48  an 

1  inch  *  2.54  an 

1  deg  lat  -  111.137  km  =  59.969  nautical  miles 

1  m/s  -  3.6  km/hr  =  1.9425  knots  =  2.2369  statute  miles/hr  =  3.2808  ft/s 
°K  -  °C  +  273.16 
°F  -  9/5  (°C)  +  32 
°Sankine  =  9/5  (°C)  +  491.69 


mw  Molecular  weight  of  water 
R*  Universal  gas  constant 
Gas  constant  for  dry  air 

g  Acceleration  due  to  gravity 
Cp  Specific  heat  of  dry  air 
pg  Standard  sea-level  pressure 

Tg  Standard  sea- level  temperature 
a  Standard  Tropospheric  lapse  rate 
k  Vcp  (approximately) 
w  Earth's  rotation  rate 


METEOROLOGICAL  CONSTANTS 
18.0  g/mole 

8.3144  X  107  ergs/moleA 

0.068557  cal/g/°K 
287.04  m2/s2/°K 

9.806  m/s2 

0.24  cal/g/°K 

1013.246  nb 
29.921  in  Hg 

288.15  °K 

0.0065  °K/m 

0.2854 

7.29X10'5/s 


K  *. 


LIST  OF  SYMBOLS  AND  ABBREVIATIONS 

a  Acceleration 

A  Altimeter  setting 

c  Wind  speed 

COL  Convective  condensation  level 
CTI  Cross  totals  index 

D  D-value  (departure  of  pressure  height  frcm  standard  atmosphere  pressure  height) 
D(t)  Displacanent  distance  as  a  function  of  time 
DD  Dewpoint  depress  ion 

dg  Spacing  (in  meters)  between  grid  points 

e  Vapor  pressure  (of  water) 

eg  Saturation  vapor  pressure 

f  Coriolis  parameter 

H  Thickness  of  a  layer  of  the  atmosphere 

H  Field  elevation 

e 

I  Grid-point  number  along  X-axis 

J  Grid-point  number  along  Y-axis 

KI  K-  index 

L  Latent  heat  of  water 

LCL  Lifting  condensation  level 

LI  Lifted  index 

LPC  Level  of  free  convection 

In  Natural  (base-e)  logarithm  (In  X  =  2.3026  log  X) 
log  Base- 10  logarithm 

P  Pressure 

PA  Pressure  altitude 

Pm  Station  pressure 

Q  Any  continuous  variable  (e.g.,  temperature,  moisture,  etc.) 
r  mixing  ratio 

rg  saturation  mixing  ratio 

Ri  Richardson  number 

RH  Relative  hunidity 

Showalter  stability  index 


SI 


reference  time 


temperature 
Dewpoint  temperature 
Wet  bulb  temperature 

Final  temperature  of  a  parcel  lifted  to  500mb 
Total  totals  index 

dx 

eastward  wind;  wind  in  X  direction  in  I,J  grid  =  ^ 
northward  wind;  wind  in  Y  direction  in  I,J  grid  =  ^ 

Vertical  totals  index 
Vertical  wind  shear 
Upward  wind 
distance  eastward 
distance  northward 

Height  of  pressure  surface  in  standard  atmosphere 
Wind  direction  (grid-relative) 

Wind  direction  (north-relative) 

Partial  differential  operator 
Finite  difference  operator 
Longitude 

Reference  longitude 
Latitude 

Potential  temperature 
Partial  potential  tenperature 
Equivalent  potential  temperature 
Pseudo-equivalent  potential  temperature 
Density  altitude 


1.  INTRODUCTION 


The  advene  of  small  computers  in  the  base  weather  station  (BWS)  has  brought  on  a  new  era  in  meteor¬ 
ological  analysis.  Among  other  uses,  conputers  permit  direct  calculation  of  many  meteorological  variables 
which  were  previously  only  estimated  from  charts,  nomograms,  or  tables.  For  some  operations,  e.g.,  com¬ 
puting  potential  temperature  (O)  for  given  temperature  (T)  and  pressure  (p) ,  there  is  only  one  equation  in 
the  literature;  it  is  familiar  to  all  meteorologists.  The  equations  used  for  some  operations  are  not  as 
well  known,  or  there  are  multiple  approximations  given  in  the  literature  (e.g.,  Tetens  formula  and  the 
Goff-Gratsch  equation  are  both  used  to  compute  saturation  vapor  pressure  (eg)  for  given  temperature) . 

This  report  lists  equations  or  algorithms  used  to  compute  armon  meteorological  and  geophysical  variables. 
It  should  serve  as  a  handy  reference  for  BWS  computer  progranmers,  promote  standardization,  and  minimize 
time-consuning  literature  researches.  Computers  can  be  applied  to  every  problem  that  involves  numerical 
calculation,  therefore,  this  list  is  not  necessarily  exhaustive.  Meteorological  algorithms  which  are 
locally  developed,  and  not  from  this  report,  should  be  submitted  to  parent  units  before  they  are  used 
operationally.  The  equations  are  given  with  little  explanation  or  interpretation,  and  users  are  cautioned 
to  consult  the  detailed  references  for  each  topic  if  more  information  is  needed.  In  most  cases,  examples 
are  given  so  local  programmers  can  validate  their  programs.  These  equations  and  algorithms  may  be  classi¬ 
fied  into  three  general  groups;  (1)  those  related  to  the  wind  fields  and  known  as  kinematic  properties; 
(2)  those  related  to  temperature  and  pressure  and  known  as  thermodynamic  properties  of  the  atmosphere; 
and  (3)  those  in  which  data  are  objectively  analyzed  to  obtain  gridded  or  contoured  data  fields.  In  this 
report,  the  meteorological  algorithms  are  provided  within  each  of  the  three  groups.  Following  these, 
relations  for  map  transformation,  chart  generation,  and  useful  geophysical  constants  are  given.  Finally, 
when  possible,  the  equations  given  here  are  taken  from  the  specifications  for  the  Automated  Weather  Dis¬ 
tribution  System  (AWDS)  so  there  will  be  minimum  discontinuity  in  data  analysis  procedures  when  AWD6 
becomes  operational. 


2.  ALGORITHMS  FDR  KINEMATIC  PROPERTIES  OF  THE  ATMOSPHERE 

2.1  Wind.  Wind  observations  are  reported  by  specifying  the  wind  direction  and  wind  speed.  For 
amputation  of  kinematic  properties  of  the  atmosphere,  it  is  necessary  to  resolve  the  wind  into  orthogonal 
components.  These  components  are  parallel  to  a  grid  system.  For  exanple,  if  a  latitude/ longitude  grid 

is  used,  the  component  parallel  to  latitude  lines  (u)  is  positive  when  the  wind  is  moving  toward  the  east, 
and  the  component  parallel  to  longitude  lines  (v)  is  positive  when  the  wind  is  moving  toward  the  north. 

The  algorithm  used  to  resolve  the  reported  direction  and  speed  to  u  and  v  components  on  a  grid  and  vice 
versa  depends  upon  the  type  of  grid  system  and  the  orientation  of  the  grid.  At  APGWC,  different  conven¬ 
tions  are  used  for  the  orientation  of  these  components  in  the  different  grid  systems.  APGWC/TH-79/003 
describes  the  convention  for  each  of  these  systems.  Units  requiring  access  to  this  technical  note  may 
request  it  through  channels  from  the  AWS  Technical  Library.  The  convention  for  wind  direction  is  to 
measure  it  clockwise  in  degrees  from  north  to  the  direction  from  which  the  wind  is  blowing. 

2.1.1  Wind  Components  from  Wind  Direction  and  Speed.  The  algorithms  for  deriving  u  and  v  components 
from  wind  direction  and  speed  are  illustrated  for  several  different  grid  systems  in  the  following: 

2. 1.1.1  Wind  Components  on  a  Latitude /Longitude  Grid.  The  u  and  v  components  are  detennined  using: 
u  =  -  c  sin  R 

v  =  -  c  cos  8 


where  8  is  the  wind  direction  and  c  is  the  wind  speed. 


2. 1.1. 2  Wind  Components  on  the  APGWC  Northern  Hemispheric  Whole-Mesh  Reference  Grid.  The  orienta¬ 
tion  of  the  u  and  v  components  and  x  and  y  axes  for  this  grid  is  given  in  Section  3. 1.2.1  and  Figure  3.6, 
AFCWC/TN- 79/003.  The  map  is  positioned  with  respect  to  the  grid  by  specifying  a  reference  longitude  A() 
which  is  parallel  to  the  horizontal  axis  (A  =  10°E) .  The  u  component  is  positive  for  a  wind  blowing 
toward  the  positive  I  direction  and  the  v  component  is  positive  for  a  wind  blowing  toward  the  negative  J 
direction.  With  this  convention,  the  u  and  v  components  of  the  wind  vector  at  a  point  are  determined 
from: 

u  =  +  c  cos  [  8  -  (  A  -  Aq)  I  2-2a 

v  =  -  c  sin  [  8  -  (  A  -  A  ) I  2-2b 

where  A  is  the  longitude  of  the  point  where  the  wind  was  observed  in  degrees  (0-360  degrees,  positive 
east  -  see  Section  2.1,  AFGWC/TN-79/003) . 

2. 1.1. 3  Wind  Components  on  the  APGWC  Southern  Hemispheric  Whole-Mesh  Reference  Grid.  The  orienta¬ 

tion  of  u  and  v  conponents  and  x  and  y  axes  for  this  grid  is  given  in  Section  3. 1.3. 2  and  Figure  3.7, 
APGWC/TN- 79/003.  The  u  component  is  positive  for  a  wind  blowing  toward  the  positive  I  direction  and  the 

v  component  is  positive  for  a  wind  blowing  toward  the  negative  J  direction.  With  this  convention,  the  u 

and  v  components  of  the  wind  at  a  point  are  determined  from: 


• 
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u  =  -  c  cos  C  3  +  (  A  -  Aq) J  2-3a 

v  =  +  c  sin  [  S  +  (  A  -  Aq) J  2-3b 

where  A  and  Aq  are  as  given  in  2.1. 1.2. 

2.1.2  Orientation  of  Winds  for  Plotting  and  Isoplething.  In  AWS,  the  convention  is  that  grid- 
relative  wind  direction  is  used  for  plotting  and  north-relative  wind  direction  is  used  for  isoplething 
(isogons).  Wind  speed  is  the  same  in  both  cases.  Given  gridded  u,v  wind  cottponents,  wind  speed  c  and 
grid-relative  wind  direction  a  are  computed  using  the  algorithm  in  2. 1.2.1.  North-relative  wind  direction 
3  is  obtained  by  sequentially  solving  the  algorithms  in  2. 1.2.1  and  2. 1.2. 2  for  the  AFGWC  grid  of  interest 
Given  north-relative  wind  direction  6,  the  algorithm  in  2. 1.2.2  can  be  used  to  solve  for  the  grid- relative 
wind  direction  i  for  the  grid  of  interest. 

2. 1.2.1  Computation  of  Wind  Speed  and  Grid-Relative  Wind  Direction  from  Gridded  Data.  All  of  the 
grids  discussed  in  Section  2.1.1  have  the  convention  that  a  clockwise  rotation  of  90  degrees  from  the 
positive  v  component  results  in  the  u  conponent.  Given  gridded  u,v  wind  conponents,  wind  speed  c  and 
grid-relative  direction  t  may  be  obtained  frcm: 


c  =  (u  +  v  )- 


=  180°  +  arctan  (u/v) 


2-4 

2-5 


where  the  sign  of  each  wind  component  is  taken  into  account  so  the  arctangent  function  returns  values  of 
:t  between  0°  and  360°  as  indicated  in  the  following  table: 


Wind  Direction  and  Speed 

Range  (  0  ^ 
a  =  180 


Sign 


180  <  i  <  270 


=  270 


270  <  »  <  360 


=  0 


0  <  .  <-  90 

<  =  90  -  0 

90  -  <  <  180  -  + 


2. 1.2.2  Relationship  Between  North-Relative  and  Grid-Relative  Wind  Directions.  North-relative  wind 
direction  /  and  grid-relative  wind  direction  a  are  related  according  to  Equations  2-6  ,  2-7  ,  and  2-8 


for  the  three  AFQWC  grids  described  above.  For  computation  of  grid- relative  wind  direction  a  fran 
observed  north-relative  direction  g,  X  is  the  longitude  of  the  reporting  station.  For  computation  of 
north-relative  wind  direction  6  from  the  grid-relative  wind  direction  a  (obtained  as  specified  in 
2. 1.2.1),  X  is  the  longitude  of  the  grid  point.  The  longitude  X  is  obtained  by  solving  Equations  3.7, 

3.8,  and  3.12,  AFQWC/TN- 79/003.  In  the  following  equations,  the  angle  solved  for  (  B  or  o,  depending  on 
the  application)  most  be  adjusted  in  miltiples  of  360°  to  keep  the  angle  between  0°  and  360° : 

AFGWC  Tropical  Grid 

B  -  a  +  180°  2-6 

AFGMC  Northern  Hemispheric  Whole-Mesh  Reference  Grid 

6  ■  a+  (90°  +  X  -  Xq)  2-7 

APGWC  Southern  Hemispheric  Whole-Mesh  Reference  Grid 

6  -  a  -  (90°  +  X  -  X  )  2-8 

2.2  Advection  and  Advected  Field. 

2.2.1  General.  Advection  is  the  process  of  transporting  an  atmospheric  property  solely  by  the 
velocity  field  of  the  atmosphere.  It  provides  an  estimate  of  the  time  rate  of  change  of  a  parameter.  The 
advected  field  is  the  field  of  values  of  the  property  after  a  given  time  interval  (during  vrtiich  advection 
is  applied). 

2.2.2  Calculations.  Mathematically,  the  horizontal  advection  of  an  atmospheric  property  Q  is: 

3Q/3t  =  -  (U3Q/3X  +  V3Q/3y)  2-9 

To  calculate  the  advection  from  values  of  U,  V  and  Q  at  grid  points  for  specified  geographical  areas, 
finite  differences  are  used  to  approximate  the  partial  derivatives  in  Equation  (2-9).  Because  of  the  use 
of  more  than  just  one  grid  type,  it  is  necessary  to  consider  the  conventions  for  positive  wind  conpcnents 
aid  positive  coordinate  axes  discussed  in  Section  2.1.  Specifically,  on  the  APCSC  polar  stereographic 
(I,  J)  grids  for  the  northern  and  southern  hemispheres  where  I  is  along  the  X  axis  and  J  is  along  the  Y 
axis,  the  instantaneous  time  rate  of  change  of  Q  at  a  point  I,  J  (advection)  is  determined  by: 

Qj.  (I,  J)  -  -  (  (U  (I.  J)  x  [Q  (I  +  1,  J)  -  Q  (I  -  1,  j)]  + 

V  (I,  J)  x  [Q  (I,  J  +  1)  -  Q  (I,  J  -  l)]}/2de  2-10 

On  the  Latitude  /Longitude  Grid,  the  instantaneous  time  rate  of  change  of  Q  at  a  point  I,  J  is  determined 
from: 

Qt  (I,  J)  -  -  {  -u  (I,  J)  x  [Q  (i  +  1,  j)  -  Q  (i  -  1,  J)]  + 

V  (I,  J)  X  [Q  (I,  j  +  1)  -  Q  (I,  J  -  1) ]}/2d  2-11 


•V-'  „• ' .  \  ‘  ■ 


• .  «  o  '  .>  .  * .  . 

-  .  V.V  V.V- . 

,  j-  '  »  v  "3*  *  •  •  *  %  ■  s*  1 


4 


vhere  Qt  (I.  J)  *  3Q/3t  is  the  time  rate  of  change  of  Q  at  point  I,  J  (advection);  u  (I,  J)  and  V  (I,  J) 
are  the  horizontal  components  of  wind  at  point  I,  J  respectively;  Q  is  the  value  of  the  atmospheric 
parameter  at  the  indicated  grid  point  (1  +  1,  J),  (I  -  1,  J),  (I,  J  +1),  and  (I,  J  -  1);  aid  dg  is  the 
grid  spacing,  a  fuiction  of  the  grid  from  which  data  are  extracted.  It  is  variable  within  the  grid.  The 
following  references  in  AFGWC/TN-79/003  allow  confutation  of  dg.  NOTE:  The  equations  in  APGWC/TN-79/003 
give  dg  in  km.  This  most  be  converted  to  the  same  distance  units  as  the  wind  speed  distance  unit.  In 
Equation  2-10,  dg  is  determined  for  the  point  (I,  J)  from  Equation  3.4,  APGWC/TN- 79/003.  The  value  of 
the  image  scale  j  in  Equation  3.4,  APCWC/TN- 79/003 ,  is  determined  from  Equation  3.5,  APGWC/TN- 79/ 003. 
The  value  of  latitude  0  in  Equation  3.5,  AP&JC/TN- 79/003,  is  determined  from  Equations  3.11,  3.7, 

and  3.8,  APGWC/TN- 79/003.  In  Equation  2-11,  dg  is  determined  for  the  point  (I,  J)  from  Equation 
3.37,  AFGWC/TN-79/003.  The  value  of  latitude  0  used  in  Equation  3.37,  APGMC/TN-79/003  is  given  by 
Equation  3.43,  APGWC/TN- 79/ 003.  The  value  of  AA  in  Equation  3.37,  APGWC/TN-79/003,  is  0.064775 
radians  and  the  value  of  AA  in  Equation  3.43,  APGWC/TN-79/003,  is  3.7113  degrees.  Suitable  approxima¬ 
tions  to  Equations  2-10  and  2-11'  cai  be  used  at  grid  boundaries  (e.g. ,  one-sided  differencing) . 

Knowing  Qt  (I,  J)  at  a  given  time  tQ,  a  predicted  value  for  Q  (advected  field)  at  a  new  time  t^  can  be 
calculated  by: 


QL  (I,  J)  =  (I,  J)  +  Qt  (I,  J)  x  (tl  -  tQ)  2-12 

where  is  the  known  value  of  Q  at  time  tQ  aid  is  the  new  value  of  Q  at  time  t^,  and  is  from 
Equation  2-10  or  2-11 .  Equation  2-12  can  be  used  to  predict  the  spatial  distribution  of  any  atmos¬ 
pheric  property  at  relatively  near  times  in  the  future  (1-12  hours) .  Equation  2-12  is  applied  in  the 
following  fashion: 


a.  The  atmospheric  property  to  be  advected  (e.g. ,  temperature)  is  specified. 


b.  The  total  time  interval  (e.g.,  current  time  plus  12  hours)  and  the  incremental  time  step  (e.g.,  1 
hour)  are  specified. 

c.  The  instantaneous  time  rate  of  change  for  the  selected  atmospheric  property  at  grid  points  is 
determined  from  Equations  2-10  or  2-11 ,  depending  on  the  grid  in  use.  The  U  and  V  components 
of  the  wind  are  assured  to  remain  invariant  in  time  for  this  computation. 


d.  Predicted  values  for  Q  are  calculated  by  Equation  2-12  for  each  time  step  until  the  total  time 
interval  is  achieved. 

2.3  Extrapolation . 

2.3.1  General.  For  short  term  forecasting  of  various  weather  features  such  as  low  and  high  pressure 
systems  and  height  and  vorticity  csiters ,  extrapolation  techniques  are  often  used  by  meteorologists.  The 
extrapolations  are  based  on  the  information  from  the  time  sequence  of  analyses  of  the  appropriate 
parameter. 

2.3.  2  Calculations.  Simple  synoptic  extrapolation  of  a  weather  feature  is  defined  by  the  fornula: 


D(t)  -  vt  +  k  at 
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where  D(t)  is  the  displacement  of  the  weather  feature  over  time,  t.  V  is  the  velocity  of  the  feature,  and 
a  is  the  acceleration  of  the  feature.  The  initial  velocity  of  the  weather  feature  is  determined  from 
confutation  of  the  distance  traveled  over  a  time  period  between  the  most  recent  known  position  at  time 
t_^  and  the  current  position  at  time  t  .  In  the  case  of  a  surface  weather  feature,  3-hour  or  6-hour 
time  intervals  of  are  most  corancnly  used.  For  rapidly  moving  features,  shorter  time  intervals  of  1  or  2 
hours  are  often  used  for  computing  initial  velocity.  For  upper  air  weather  features,  the  initial  velocity 
is  determined  over  a  12-hour  interval  between  the  currant  data  observation  time  and  12  hours  earlier. 

Thus,  the  time  interval  for  confuting  initial  velocity  is  variable  and  should  be  an  input  parameter.  To 
confute  acceleration,  the  position  before  t_^  (t_2^  *-s  a^so  required.  The  technique  for  determining 
velocity  and  acceleration  can  use  a  cartesian  coordinate  system.  The  current  position,  XQ  at  time  t  ,  and 
past  two  positions,  X_^  at  time  t_^  and  X_2  at  time  t_2  are  determined  in  x,  y  coordinates.  If  only  two 
positions  of  a  feature  are  available,  the  method  may  be  applied  by  assuming  that  the  acceleration  is  zero. 
To  apply  the  extrapolation  technique,  compute  displacements  along  the  *  and  y  axes  separately.  Equation 
2-14  can  be  used  to  confute  the  extrapolated  position  along  the  X-axis  and  Equaticn  2-15  can  be  used 
to  compute  the  extrapolated  position  for  the  Y-axis  (X+^  and  Vl  respectively) : 

x+i  -  xo  +  [(x0  -  x_1/(tQ  -  t^)  i  <t+1  -  tQ)  + 

C(X0  -  x_1)/(tQ  -  t_x>  -  (x_j_  -  x_2)/(t.1  -  t_2)jr(t+1  -  to)2/(to  -  t_2)j  2-14 

Y+1  =  Yo  +  [<Yo  -  Y-l)/(to  -  C-l>]  <c+l  '  Co>  + 

[(Yo  -  Y-l)/(to  -  C-l>  -  <Y-1  -  Y-2)/(t-l  -  t-2>  l[(t+l  -  -  C-2>  1  2-15 

2.4  Vertical  Wind  Shear.  The  vertical  wind  shear  (VWS)  is  the  difference  between  wind  vectors 
between  two  levels  in  the  atmosphere.  The  magnitude  of  the  VWS  is  usually  all  that  is  required  for 
analysis  and  is  confuted  from  reported  or  interpolated  u  and  v  wind  conponents  using  the  following  rela¬ 
tionship  : 

VWS  =  [(U2  -  Uj)2  +  (V2  -  Vj)^]^  2-16 

where  the  subscripts  1  and  2  refer  to  the  lower  and  upper  atmospheric  levels,  respectively.  When  reported 
wind  speed  (c)  and  direction  (6)  are  available,  the  raagiitude  of  the  VWS  is  computed  using  the  following 
relationship: 

2 

VWS  =  [(c2  sin  82  ~  sin  8^)  + 

(C2  cos  82  -  cos  8j)2]^  2-17 

where  the  subscripts  1  and  2  refer  to  the  lower  and  upper  atmospheric  levels,  respectively. 

2.5  Streamline  Analysis.  Streamline  analyses  are  generated  using  wind  direction  and  wind  speed  data 
as  input  within  selected  regions .  Following  a  procedure  described  by  Whittaker  (1977) ,  in  Monthly  Weather 
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Review  (Volune  105,  pp.  786-788),  the  slopes  of  the  streamlines  can  be  determined  from  observed  wind 
directions .  The  display  can  consist  of  a  series  of  lines  tangent  to  die  reported  wind  directions  (within 
an  instantaieous  wind  flow  pattern).  Ihits  requiring  access  to  Whittaker's  article  should  forward 
requests  to  the  AWS  Technical  Library  through  their  parent  unit. 


3.  THERM3DTOAMIC  VARIABLES  AND  THEIR  ALGORITHMS 

3.1  General .  Observations  of  the  atmosphere  above  the  ground  are  taken  at  selected  station 
locations  over  the  globe  four  times  each  day  at  0000,  0600,  1200,  and  1800  GMT  and  twice  each  day  at  0000 
and  1200  CMT  at  other  selected  stations.  These  observations  are  input  to  the  derivation  of  thermodynamic 
variables  that  are  used  as  diagnostic  and  forecasting  tools.  Temperature,  pressure  (or  height) ,  and 
humidity  are  measured  at  the  surface  and  aloft.  Using  the  First  Law  of  Thermodynamics,  the  equation  of 
state,  Dalton's  Law,  and  other  fundamental  relations,  a  variety  of  thermodynamic  variables  can  be  derived. 
Meteorologists  in  the  BWS  have  usually  made  these  analyses  by  plotting  the  observed  data  on  charts,  such 
as  the  Skew-T,  Log-P  diagran,  and  by  estimating  other  variables  from  lines  on  the  chart.  Of  course,  the 
derived  quantities  can  be  ccnputed  directly  frcrn  the  observations  if  desired.  Equations  and  algorithms 
for  some  oomnon  thermodynamics  variables  are  given  below.  Uhless  stated  otherwise,  the  equations  listed 
are  those  given  in  the  AWDS  specifications.  For  application  of  these  thermodynamic  quantities,  fore¬ 
casters  should  consult  other  appropriate  reference  documents  (e.g.,  AWS/TR- 79/006  discusses  general  theory 
and  use  of  the  Skew-T,  Log-P  diagram;  AWS/TR-80/001  discusses  aircraft  icing  forecasts,  etc.). 

3.2  Glossary  of  Terms.  The  following  is  a  glossary  of  terms  that  appear  in  this  section.  Part  A 
consists  of  terms  that  are  derived  from  observed  parameters.  Part  B  lists  the  observed  paraneters  and 
Part  C  contains  other  parameters,  some  of  which  are  ccnputed  from  the  observed  parameters  that  are  input 

to  the  computation  for  the  derived  parameters,  and  others  which  are  constants.  The  equations  for  ccnputing 
the  parameters  in  Part  C  are  given  once  here,  rather  than  being  repeated  several  times  throughout  the  text 
of  this  report.  Part  D  contains  useful  conversion  factors  for  the  computations . 

Part  A.  Derived  Parameters: 

0,  Potential  temperature  (°K) 

0,  Mean  potential  temperature  of  a  layer  (°K) 

H,  Thickness  of  a  layer  of  the  atmosphere  (m) 

D,  D- values  (m)  [departure  from  the  standard  atmospheric  height  of  a  constant  pressure 

surface] 

r,  Mixing  ratio  (gjn/kg) 

?,  Mean  mixing  ratio  of  a  layer  (gn/kg) 

r  ,  Saturation  mixing  ratio  (gn/kg) 

RH,  Relative  humidity  (%) 

LCL,  Lifting  condensation  level 

Temperature  at  the  LCL  (°K) 


Pressure  of  the  LCL  (mb) 


■  1^.  !  ■,■  I  I  I T i  iimm 
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Height  of  the  LCL  (m) 

CCL,  Convective  condensation  level  (mb) ] 

LFC,  Level  of  free  convection  [Py^  (mb)] 

SI,  Showalter  index  (°K) 

LI,  Lifted  index  (°K) 

KI,  K- index  (°K) 

VTI,  Vertical  totals  index  (°K) 

CTI  Cross  totals  index  (°K) 

TTI,  Total  totals  index  (°K) 

Part  B.  Observed  Parameters  Used  for  Computation  of  Derived  Parameters: 

P,  Atmospheric  pressure  (mb) 

Z,  Height  of  a  constant  pressure  surface  (m) 

T,  Dry  bulb  tenperature  (°C) 

DD,  Dewpoint  depression  (°C) 

A,  Altimeter  setting  (in,  Hg.) 

Part  C.  Constants  and  other  Input  Parameters  to  Required  Derived  Parameters  (Part  A)  Computed  from 
Observed  Parameters : 

C.l  Constants  and  Other  Variables: 

nw,  Molecular  weight  of  water  vapor  -18.0  gn/mole 

R*,  Universal  gas  constant  «  8 . 3144  «  10^  ergs/nDle/°K 

R^,  Gas  constant  for  dry  air,  -  0.068557  cal/gjn/°K 

-  2.8704  «  102  ra2/sec2/°K 

2 

g,  Acceleration  of  gravity,  g  -  9.806  m/sec 

Cp,  Specific  heat  of  dry  air  at  constant  pressure  «  .24  cal/gjn/°K 

Zg.  Height  of  a  constant  pressure  surface  in  the  standard  atmosphere  (m) 
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P  ,  Standard  sea  level  pressure  -  1013.246  nb  or  29.921  inches  Hr.  or 
101.325  kilopascals  (kpa) 

T  ,  Standard  sea  level  tenperature  =  288.16°K 

P  .  Station  pressure  (mb) 

m 

PA,  Pressure  altitude  (m) 

H’  Density  altitude  (m) 

a,  The  standard  atmospheric  lapse  rate  below  the  isothermal  layer,  a  =  0.0065°  C/m 

k.  Ratio  of  the  gas  constant  for  dry  air  to  the  specific  heat  capacity  of  dry  air  at 

constant  pressure  =  0.2854 

H  ,  Field  elevation  (m) 

Tenperatures  are  reported  in  °C,  but  except  where  otherwise  indicated,  temperature  in  the  equations 
in  the  following  sections  use  °K,  where  T(°K)  =  T(°C)  +  273.16. 

C.2  Other  Input  Parameters: 

(a)  T^,  Dewpoint  temperature  (°C),  upper  air 

Td  =  T  -  DD  (°C)  C-2(a) 

(b)  e,  Vapor  pressure^  ^  (mb) 
when  T  >  273.16  °K,  then 

e  -  10**  [23.832241  -  5.02808  *  Log  <Tj) 

-  1.3816  x  (10**(-7))  x  (10**(11. 334  -  0.0303998  x  Tj)) 

+  8.1328  x  (10**(-3))  x  (10**(3. 49149  -  1302 . 8844/Td) ) 

-  2949.076/Tj]  C-2(b) 

and  when  T  <  273 . 16  °K,  then 

(t)  The  relation  given  is  the  Goff-Gratch  formila,  which  is  correct  for  all  values  of  Tj. 

However,  as  the  Goff-Gratch  formula  is  often  ember  some,  approximate  formulae  are  sometimes  used  (a 

detailed  discussion  is  given  by  David  Bolton,  Monthly  Weather  Review,  Volume  108,  pp.  1046-1053).  One 

formula  with  accuracy  suitable  for  most  purposes  is  Tetens'  formula:  ^  a^d  ^ 

e  -  6. limb  x  10**  (Td  +  b) 

where  a  -  7.5,  b  -  237.3  for  Td  >  0°C,  and  a  -  9.5,  b-  265.5  for  Td  <  0°C. 


C-2(b) (1) 


e  -  10**  [3.56654  x  Log  (Td)  -  0.0032098  x  Td 


-  2484.956/Td  +  2.0702294] 


C-2(c) 


where  Td  is  dewpoint  temperature  in  °K, 


**  means  "raise  to  the  exponent  of," 
Log  means  "take  to  base-10  logarithm" 
Check  values  for  these  equations  are: 


Equations  C-2(b)  or  C-2(c) 

Equation  C-2  (b)(1) 

50°C 

e  =  123.40nb 

e  -  123.39mb 

20°C 

23.373 

23.389 

-10°C 

2.597 

2.596 

-40°C 

0.1283 

0.1262 

(c)  eg,  Saturation  vapor  pressure 


To  compute  e  ,  use  Equations  C-2(b)  and  C-2(c),  C-2(b)(l),  replacing  the  dewpoint  with 
s 


with  the  temperature . 

(d)  L,  Latent  heat  of  water  vapor  (cal/gpj) 

L  -  597.3  -  0.564  (T  -  273.16)  when  T  >  273.16°k 

L  -  597.3  -  0.574  (T  -  273.16)  when  T  <  273.16°K 


C-2(d) 

C-2(e) 


(e)  ege>  Pseudo-equivalent  potential  tenperature  (°K) 


9se  *  9d  (Lr/cpT) 


C-2(f) 


where  "exp"  means  raise  the  natural  base  e  to  the  exponent  in  the  argunent  that 
follows,  and  r  is  the  mixing  ratio  in  grams  per  gram  (given  in  g/Kg  by  3-2)  and 
ed  is  given  by  C-2(g) . 

(f)  9d,  Partial  potential  tenperature  (°K) 

9d  -  T  [(1000)/(P  -  e)]k  C-2(g) 
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(g)  T,n  (parcel),  Final  estimated  tenperature  of  a  parcel  lifted  to  500  mb. 


Part  D.  Conversion  Factors: 


1  mb  -  1000  dynes/cnf 
1  in.  (Hg)  -  33.86389  nb 


1  dyne  »  1  gn  cm/sec 


1  foot  »  . 3048  m 


1  erg  «  1  dyne  an  »  2.3892  *  10_o  cal. 

1  calorie  *■  4.186  x  10^  ergs 
=  4.186  joules 

3.3  Potential  Temperature,  Mixing  Ratio,  Saturation  Mixing  Ratio,  Relative  Humidity,  and  Layer 


Average  of  Quantity  Q. 

3.3.1  Potential  Temperature.  The  potential  temperature  (0)  is  the  temperature  a  parcel  of  air  would 
have  if  brought  dry  adiabatically  to  a  standard  pressure  of  1000  mb. 

0  -  T  [(1000/P)k] 

A  check  value  is:  if  P  *  850  and  T  =  20°C  then  0  *  307.1^  3-1 


3.3.2  Mixing  Ratio.  The  mixing  ratio,  r  is  defined  by: 


r  =  1000  x  . 622e/(P-e) 


vbere  e  is  given  by  Equation  C-2(b)  or  C-2(b)(l)  or  C-2(c).  A  check  value  is: 
if  Tj  =  20°C  and  P  =  990  mb  then  r  =  15-04g/kg. 

3.3.3  Saturation  Mixing  Ratio.  The  saturation  mixing  ratio,  rg  is  defined  by; 

rg  -  1000  x  .622eg/(P-eg) 

tdiere  eg  is  given  by  Equation  C-2(b)  or  C-2(b)(l)  or  C-2(c).  A  check  value  is: 

if  P  =  980  nb  and  e„=  20  nb  then  r  =  12 . 958 . 

s  s 

3.3.4  Relative  Humidity.  The  relative  humidity,  RH  at  any  level  is  computed  from: 

RH(%)  -  r  (100) 


v-  v  v  /-  -  z-y  yy  y-y  yy-  •  -w  ■  -  .  .  .  .  .  -  - .  • .  * .  - 

A  1  mm jCAm—L  hi  .  iimldl.eilii.i  <  ■■  ■  ■  ■«  W  ■  il««i  ill  <  n,.-.  l 


.*  «• 


; 


where  r  and  rg  are  determined  from  Equations  3-2  and  3-3  ,  above.  A  check  value  is: 
if  T  =  20°C  and  Tj  =  10°C  and  P  =  990  mb  then  RH  =  51.9%. 

3.3.5  Layer  Average  of  Quantity  Q.  The  average  of  a  quantity  Q  is  calculated  from  Equation  3-5: 


n-1 


q  (Pm,P'>  =  [1/ (In  Pm  -  In  P1)]  X  [(1/2)  X  (Q.  +  Q.+1)  X 


(In  P.  -  In  P.+1)J  +  [(1/2)  x  (^  +  Q')  *  (In  Pn  -  In  P1)  ]  } 
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where 


Q  =  average  of  Q  over  the  pressure  layer  from  to  P' 
Pm  =  station  pressure  (pressure  at  bottom  of  the  layer) 
=  value  of  Q  at  Pm 
P'  =  pressure  at  top  of  layer 


Q^,  =  values  of  Q  and  pressure  respectively  at  level  i 

P  =  highest  level  i  such  that  Pm  >  P  >  P' 
n  **  m  —  n 


<!'  •  %  *  <Vl  -  V  *  «”  Pn  -  P„  -  W 

where  "In"  means  "take  the  natural  logarithm." 

The  first  term  of  the  two  terms  in  curly  brackets  in  Equation  3-5  is  zero  if  n  =  m  (i.e. ,  if  Pm  =  Pn  >  P') 

3.4  Equivalent  Potential  Temperature .  The  equivalent  potential  temperature  (0g)  takes  into  account 
the  possible  gain  (loss)  of  sensible  heat  that  takes  place  during  condensation  (evaporation)  and  is 
defined  by  the  dry  adiabatic  reduction  of  the  equivalent  temperature  to  a  standard  pressure  of  1000  mb. 

9g  is  computed  from: 


0g  =  T  [(1000/P)  k<1-°-28r)]X 

exp[(3376/TLCL  -  2.54)r(l  +  0.81r)],  3-6 


vbere  K  =  0.2854,  P  is  the  pressure  in  nb,  T^^  is  the  temperature  in  °K  at  the  lifting  condensation 
level  given  by  Equation  3-9  and  r  is  the  mixing  ratio  in  grans  per  gran.  A  check  value  is: 


if  T 
LCL 


=»  20°C,  Td  =  10°C,  and  P  =  995  mb,  then  eg 

»  280.97°K,  0  -  316.14°K. 
e 


12.272  mb,  r  =  7.767  x  10"3  g/g, 
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3.5  Thickness  and  D- Values . 


3.5.1  Thickness.  Thickness,  H  is  computed  from: 

H  =  Z2  -  Zx  3-7 

Z^  is  the  height  of  the  lower  (in  height)  constant  pressure  surface  and  Z2  is  the  height  of  the  higher 
constant  pressure  surface.  Examples  of  constant  pressure  surface  layers  for  which  thickness  fields  are 
computed  and  their  thickness  in  the  standard  atmosphere  are: 

Layer  Thickness  in  the  Standard  Atmosphere 


a. 

1000  -  500  mb 

5463m 

b. 

1000  -  700  mb 

2901 

c. 

1000  -  850  mb 

1346 

d. 

850  -  500  mb 

4117 

e. 

850  -  700  mb 

1555 

3.5. 

.2  D- Values.  The  D-value  is  computed  from: 

D  =  Z  -  Zg  3-8 

Zg  is  the  height  of  a  pressure  surface  in  the  standard  atmosphere  and  Z  is  the  reported  height  of  that 
pressure  surface.  Table  3-1  gives  the  standard  height  (Z  )  for  each  standard  constant  pressure  surface 


Table  3-1 


Standard  Atmosphere  Heights  of  Constant  Pressure  Surfaces 


ESS 'JRE 

HEIGHT 

1  (nb) 

Zs  (m) 

1000 

111 

950 

540 

900 

988 

850 

1,457 

800 

1,949 

750 

2,466 

700 

3,012 

650 

3,591 

600 

4,206 

550 

4,865 

500 

5,574 

450 

6,344 

400 

7,185 

350 

8,117 

300 

9,164 

250 

10,363 

200 

11,784 

150 

13,608 

100 

16,180 

70 

18,442 

50 

20,576 

30 

23,849 

D-values  can  be  conputed  for  any  or  all  of  the  constant  pressure  surfaces. 

3.6  Lifting  Condensation  Level.  The  lifting  condensation  level,  (LCL) ,  is  the  height  at  which  a 
parcel  of  air  would  became  saturated  when  lifted  dry  adiabatically.  The  temperature  of  the  LCL,  is 

calculated  from: 

TLCL  =  Td  "  (°-212  +  0.00157lTj  -  0.000436  T  )  fT  -  Td)  +  273.16  3-9 

where  Tj  is  the  surface  dewpoint  temperature  and  Tfl  is  the  surface  temperature ,  both  in  °C.  A  check  value 
is: 


if  Tq  =  20°C  and  Td  =  10°C  then  TLa  =  280.97°K. 


The  pressure  of  the  LCL  is  given  by: 
PLCL  ’  Po 
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where  PQ  is  the  surface  pressure  and  TQ  is  the  surface  temperature  in  °K.  A  check  value  is : 

if  tQ  =  20°C,  PQ  =  980  mb,  =  281°K  then  P^  =  344.3  nb 
The  height  of  the  LCL  is  given  by: 

-  (PT/g)  In  (P^P^)  3-11 

T  is  computed  from  Equation  (3-5)  with  T  (P  ,  PTrT)  replacing  Q  (P  ,  P')  in  Equation  3-5- 

o  LaJj  m 

3.7  Convective  Condensation  Level.  The  convective  condensation  level,  CCL,  is  defined  as  the  level 
at  which  air  that  is  heated  near  the  earth' s  surface  will  rise  ndiabatically  until  it  becomes 
saturated  (condensation  begins).  To  determine  the  CCL: 

(a)  Calculate  the  average  mixing  ratio  (r)  in  the  layer  from  the  surface  to  100  trb  above  the 
surface  from  Equations  3-2  and  3-5  replacing  Q  (P^,  P')  in  Equation  3-5  with  r 
(PQ,  PQ-100).  In  the  following,  the  average  mixing  ratio  in  the  lowest  100  mb  will  be 
denoted  by  r. 

(b)  Calculate  the  saturation  mixing  ratio  rg(i)  from  Equation  3-3  at  each  pressure  level 
reported,  until  rg(i)  at  a  given  pressure  level  is  less  than  or  equal  to  r. 

(c)  If 

(1)  rg(i)  =  ?,  then  Pca  =  P(i) 

(2)  rg(i)  <  r,  caipute  PCCL  by  interpolating  between  that  level  and  the  Tcal  1  wer  level: 

Pca  =  -  ^(P(i)  '  P(i-l»/(r8(i-l)  -  rg(i)'  '  +  P(i-l)  3-12 

3.8  Level  of  Free  Convection.  The  level  of  free  convection,  (LPC) ,  is  defined  as  the  height  at 
which  a  parcel  of  air  lifted  dry  adiabatically  until  saturated  and  moist  adiabatically  thereafter  would 
become  warmer  (less  dense)  than  its  environment.  The  procedure  for  calculating  the  pressure  Py^,  of  the 
LPC  includes  an  approximate  numerical  procedure  which  is  given  at  the  end  of  this  section.  The  actual 
calculation  of  the  pressure  of  the  LPC  proceeds  as  follows: 

(a)  Calculate  the  temperature  and  pressure  of  the  LCL  from  the  algorithms  in  Section  3.5. 

(b)  Calculate  the  pseudo-equivalent  potential  tarperature  Gge  at  the  LCL  from  Equations 

C-2(f)  and  C-2(g)  ,  replacing  T  with  Ty^  and  P  with  PLCL  NOTE:  Above  the  LCL,  the 
parcel  is  saturated;  therefore,  its  temperature  and  dewpoint  temperature  are  equal.  Pro¬ 
ceed  to  step  (c)  with  P(i)  taken  as  the  lowest  height  level  where  P(i)  <  Py^. 


(c)  Perform  the  approximate  numerical  procedure  with  e  -  0.05  to  solve  for  T^(i) . 

(d)  Confute  the  difference  A  T(i)  between  the  reported  tenperature  Tg(i)  and  the  ccnputed 
tenperature  Tp(i) : 


A  T(i)  -  T  (i)  -  T  (i) 
a  p 


If  A  T(i)  =  0,  then  take  Tp(i)  as  the  tenperature  T^^  of  the  LPC  and  proceed  to  step  (h). 

(e)  If  at  step  (d) ,  A  T(i)  <  0,  proceed  to  step  (h)  after  computing  with  the  following 
equation: 


LFC 


[  AT(i)  T  (i-1)  -  AT(i-l)T  (i)]/[  AT(i)  -  AT(i-l)  ] 
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(f)  If  at  step  (d),  AT(i)  >  0,  increment  the  pressure  level  counter  i  (go  to  the  next  higher 
height  level).  If  i  <  n,  where  n  is  the  nunber  of  reported  levels  in  the  sounding,  return 
to  step  (c) . 

(g)  If  at  step  (f),  i  >  n,  terminate  the  calculations,  and  indicate  that  the  LFC  does  not  exist. 

(h)  Perform  the  approximate  numerical  procedure  with  e  =  0.1  to  solve  for  Pjj^,  replacing 
Tp(i),  P(i) ,  T,  T',  and  AT  with  T^p^,  Pjjq.  P‘ ,  and  AP,  respectively. 

The  approximate  numerical  procedure  uses  a  quantity  E: 


E  =  Tp(i)  [  (1000)  /  (P(i)  -  es)]K  expCLr,/CTp(i)  ]  -  0ge 
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The  objective  is  to  have  E  be  as  close  to  zero  as  possible  to  obtain  the  closest  approximation  for  the 
desired  parameter.  The  following  procedure  is  carried  out: 

(a)  Estimate  a  value  of  T  and  an  arbitrary  parameter  value  increment  AT. 

(b)  Calculate  E  from  Equation  3-14  and  test  for  |E|  <  e,  where  "|  |"  means  "take  the 
absolute  value  of."  If  E  passes  this  test,  then  T  is  taken  as  tine  desired  parameter 
value.  If  not, 

(c)  Form  T'  =  T  +  AT  and  an  associated  E' .  Again,  test  for  | E ' |  <  e.  If  E'  passes  the  test, 
then  T'  is  taken  as  the  desired  parameter  value.  If  not, 


(d)  Compare  the  signs  of  E  and  E' .  If  they  are  different  (i.e.,  the  desired  zero  crossing  of 
E  is  between  T  and  T'),  then  divide  AT  by  2  and  return  to  step  (c). 


(e)  If  at  step  (d) ,  the  signs  of  E  and  E'  were  the  same,  compare  the  absolute  values  of  E  and 
E' . 

(f)  If  | E’  |  <  [E| ,  set  T  equal  to  the  previous  value  of  T' ,  and  E  equal  to  the  previous  value 
of  E' ,  and  return  to  step  (c) . 

(g)  If  [E'|  >  |E|  change  the  sign  of  AT  and  return  to  step  (c). 


3.3.1  Dewpoint  Temperature  from  Wet  Bulb  Tenperature.  In  same  instances  the  wet  bulb  tenperature 
(T) ,  such  as  measured  with  a  sling  psychrometer,  may  be  available  instead  of  the  dewpoint  temperature 


(Tj).  The  following  procedure  may  be  used  to  compute  T^.- 

(a)  Compute  0Q  =  TQ(1000/po)  '2^  where  Tq,  Pq  are  surface  temperature  aid  pressure. 

(b)  As  in  Section  3.8,  confute  0ge  from  Equations  C-2(f)  and  C-2(g)  using  Twand  T 

(c)  Follow  the  procedure  described  in  Section  3.8  to  compute  ^  and  P.  at  levels  above  the 

surface.  At  each  level  i,  conpute  0.  using  T.  and  P.  and  Equation  3-1.  If  ,0  -  0  I 

1  11  1  o 1 

0  0-  «,  step  (d)  belcw.  If  0^  <  0^,  procetu  ■  ne  next  higner  level.  If  0^  >  0Q 

reduce  the  size  by  vdiich  P  is  changed  by  half  and  recompute  T. ,  P. 

(d)  Td  (°C)  is  given  by 

T.  =  318.27  -  0.002007  T  -  1.212  +  SQ 

Q  O 


where 


SQ  =  (1.212  -  0.002007  Tq)2  - 

0.006284  [273.16  -  T.  -  T  (0.121  -  0.000436  Tq)]  1,2 

and  T  is  °C  and  T.  is  °K. 
o  1 

3.9  Stability  Indices.  Stability  indices  are  used  to  provide  a  quantitative  basis  for  the  predic¬ 
tion  of  showers,  thunderstorms,  and  severe  local  storms.  Calculations  for  several  of  the  stability 
indices  are  described  belcw. 

3.9.1  Showalter  Stability  Index.  Showalter  developed  a  stability  index  (SI)  in  1953  that  is  widely 

used.  It  is  the  difference  between  the  free  air  temperature  T^r,  at  the  500  mb  level  and  the  temperature 
A  3U 

T50  (parcel)  of  a  parcel  of  air  that  is  lifted  dry  adiabatically  from  850  nib  to  its  proper  LCL,  then 
moist  adiabatically  to  500  mb.  To  determine  the  SSI: 

(a)  Determine  the  temperature  T^^  t^le  Pressure  plcL,  t*le  above  850  mb  from 
Equations  3-9  aid  3-10  ,  replacing  TQ  and  T^  with  T85  and  Td85  (the  temperature  and 
dewpoint  temperature  at  pressure  Pg5  at  850  mb) ,  and  setting  PQ  -  Pg^ . 

(b)  Determine  the  pseudo-equivalent  potential  tenperature  0ge  at  the  LCL  (above  850  mb)  from 
Equations  C-2(f)  and  C-2(g)  ,  replacing  T  with  T^^  a11'1  p  with  Pj^.  NOTE:  Above  the 
LCL,  the  parcel  is  saturated;  therefore,  its  tenperature  and  dewpoint  taiperature  are  equal. 

(c)  Perform  the  approximate  numerical  procedure  given  in  Section  3.8  with  -  0.05  to  solve  for 

A  A 

T50  (parcel),  replacing  T^(i)  and  P(i)  with  T^g  (parcel)  and  P^g  *>  500  mb,  respectively. 

(d)  Calculate  the  SSI  from  the  following  relationship: 


A 

SSI  -  Tcn  -  Tcn  (parcel) 
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3.9.2  Lifted  Index.  There  are  times  when  the  850  mb  parameters  are  not  representative  of  the  tem¬ 
perature  and  molstjtrp  omHitions  in  the  lower  atmosphere,  thereby  causing  the  SSI  to  be  unrepresentative 
of  the  stability.  To  overcame  this  problem,  the  Lifted  Index  (LI)  has  been  developed.  To  determine  the 
LI: 


(a)  Determine  the  temperature  T^^  and  the  pressure  PL(^  at  a  modified  LCL  from  Equations 

3-9  and  3-10  ,  replacing  P  with  Pq  -  50  trb,  and  replacing  Tq  and  T^  with  T  and  f ,, 
the  average  temperature  and  dewpoint  temperature  frcm  the  surface  to  100  mb  above  the 
surface.  The  quantities  T  and  are  computed  from  Equation  3-5  ,  replacing 

Q(Pm,  P')  with  T(Pq,  P  -  100)  and  Td(PQ,  PQ  -  100). 

(b)  Determine  the  pseudo-equivalent  tanperature  6se  at  the  modified  LCL  from  Equations 
(C-2(f))  and  (C-2(g)),  replacing  T  with  T^^  and  P  with  PTrT .  NOTE:  Above  the  LCL,  the 
parcel  is  saturated;  therefore,  its  temperature  and  dewpoint  tanperature  are  equal. 

(c)  Perform  the  approximate  numerical  procedure  given  in  Section  3.8  with  e  =  0.05  to  solve  for 

A  A 

T50  (parcel),  replacing  T^(i)  and  P(i)  with  T^q  (parcel)  and  P^q  =  500  mb,  respectively. 

(d)  Calculate  the  LI  from  the  following  relationship: 

A 

LI  =  T50  -  T50  (parcel)  3-16 

3.9.3  K- Index.  The  equation  for  calculating  the  K- index  is: 


KI 


T85  +  Td85  "  T70  +  Td70  "  T50 
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where 


Tgs  =  850  mb  temperature  (°C) . 

Td85  =  ®^0  dewpoint  tanperature  (°C) . 

T?0  =  700  mb  temperature  (°C) . 

T^q  =  700  mb  dewpoint  temperature  (°C). 

T^q  =  500  mb  temperature  (° C) . 

3.9.4  Vertical  Totals  Index.  The  equation  for  calculating  the  Vertical  Totals  Index  (VTI)  is: 

vn  =  t85  -  t50  3-18 

where 


T 


85  " 


T 


cn  = 


850  mb  temperature  (°C) . 
500  mb  temperature  (°C) . 


3.9.5  Cross  Totals  Index.  The  equation  for  calculating  the  Cross  Totals  Index  (CTI)  is: 

“  Td85  -  T50  3'19 

where 

^d85  *  “b  dewpoint  temperature  (°C) . 

Tjg  =  500  mb  tenperature  (°C) . 

3.9.6  Total  Totals  Index.  The  equation  for  calculating  the  Total  Totals  Index  (1TI)  is: 

TTI  =  VTI  +  CTI  3-20 

where 

VTI  =  vertical  totals  index. 

CTI  -  cross  totals  index. 

3.9.7  Richardson  Nunber .  The  Richardson  matter  (Ri)  is  a  nondimensional  ratio  of  static  stability 
and  wind  shear.  It  is  usually  computed  across  vertical  layers  less  than  about  5000  feet  in  depth.  Values 
less  than  0.25  indicate  turbulence  aid  values  0.25  to  1.0  indicate  a  chance  of  turbulence.  The  equation 
for  Ri  is: 

Ri  =  g(02  -  01)  H/(VWS2)0 

2 

where  g  =  9.806  m/s  ,  02  and  0^  are  the  potential  temperatures  at  the  top  and  bottom  of  the  layer 
respectively,  H  is  the  thickness  from  Equation  3-7  ,  ¥  =  (0^  +  02)/2,  and  VWS  is  the  vertical  wind  shear 
from  Equations  2-16  or  2-17.  A  check  value  is:  if  0j^  =  347°K,  02  =  353°K,  H  -  150Gta,  and  VWS  =  10 
m/s,  then  Ri  =  2.521. 

3.10  Density  Altitude.  The  altitude  above  sea  level  characterized  by  a  known  air  density  is  called 
the  density  altitude.  It  is  computed  from  Equation  3-21  and  this  equation  is  valid  up  to  a  height  of 
35,332  feet: 

f ^  =  145,366  [1  -  (17.326  P/y235)  3^21 

where 

*  density  altitude  in  feet 

Tv  -  -0.288  +  (9/5)  T  (1  +  1.60779r)/(l  +  r)  3-22 

where 

Tv  is  (°  Raikine) ,  T  is  (°K) ,  r  is  (gjn/gn) ,  and 

P  “  Station  barometric  pressure,  in  inches  of  mercury: 


'3-23 


Pm  “  tAVn  - 

where  A  -  altimeter  setting  (inches  of  mercury) , 


n  -  5.2561, 

Hg  =*  field  elevation  (m), 

and  a,  Pg  and  Tg  are  listed  in  Section  3.2.C.I.  A  check  value  is:  if  Hg  -  105m, 

A  =  30.02  im. ,  T  =  283. 2°K,  r  =  .006  g/g,  then  Ty  -  511.32°R  and  Pm  -  29.648  and  ^  -  141.43  feet. 

3.11  Pressure  Altitude.  The  pressure  altitude  (PA)  is  the  height  above  sea  level  in  the  standard 
atmosphere  at  vhich  the  altimeter  setting  (A)  exists.  It  is  computed  from  Equation  3-2y-  : 

PA  =  (Tg/a)  [1  -  (A/Pg)1/n]  +He  3-24 

where 

A  =  Altimeter  setting  (reported) 

PA  =  Pressure  altitude  for  the  Altimeter  setting 
n  =  5.2561 

A  check  value  is:  if  A  =  30.02  and  Hg  =  105m,  then  PA  =  77.13  m. 

3.12  Sea-Level  Pressure.  The  procedures  used  within  AWS  to  reduce  station  pressure  to  sea- level 
pressure  (SLP)  are  described  in  IMH  8(b) .  As  the  required  procedures  involve  using  a  12-hour  "mean" 
tenperature  and  the  pressure  reduction  (hand-held)  computer  (WBAN  54-7-8  or  CP-402 /UM) ,  together  with  a 
table  of  pressure  reduction  ratios,  there  seems  to  be  little  advantage  gained  at  this  time  by  using  a 
minicomputer  for  this  process. 

4.  SPATIAL  ANALYSIS  ALGORITHMS 


4. 1  General.  Currently  analyses  of  meteorological  parameters  at  base  weather  stations  are  either 
done  manually  or  they  are  received  via  the  Air  Force  Digital  Graphics  System  (AFDIGS)  from  AFGMC.  Local 
units  may  desire  to  have  the  capability  for  automated  objective  analyses  residing  at  the  BW5  for  both 
those  currently  produced  manually  and  those  analyses  received  from  APGWC.  Before  they  can  perform  the 
analyses,  they  mist  first  transform  observations  at  irregularly  spaced  points  into  data  at  points  on  a 
uniformly  spaced  grid;  i.e.,  generate  horizontal  grid  data.  The  automated  analyses  that  may  be  generated 
at  the  BWS  include  the  surface  and  upper  air  local  Area  Work  Chart  (LAWS) ,  and  cross  section  analyses  for 
various  parameters.  The  LAMC  analyses  are  in  the  horizontal  (x,  y)  plane  and  the  cross  sections  are  in 
the  vertical  (log  P)  plane. 


Horizontal  Objective  Analysis.  The  "nearest  neighbor"  technique  described  in  Section  4.2.1 
(Western  Region  Technical  Attachment  No.  82-14)  can  be  used  to  develop  grid  values  from  station  data. 
Grid  points  are  assigned  values  of  the  closest  observation  to  the  grid  point. 
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4.2.1  The  Nearest  Neighbor  Analysis  Technique.  The  generation  of  grid  point  values  using  the 
nearest  neighbor  technique  is  basically  a  three  step  procedure: 

(a)  Step  1.  Assisi  a  value  at  a  grid  point  by  setting  it  equal  to  the  value  of  the 
observation  closest  to  the  grid  point.  Each  observation  can  be  assigied  to  only 
one  grid  point.  By  using  observations  only  once,  the  distance  from  the  observation 
to  the  assigned  grid  point  is  less  than  one-half  of  the  grid  spacing.  If  there  is 
more  than  one  observation  within  one-half  grid  interval  of  any  grid  point,  then  the 
observation  farthest  from  the  grid  point  may  be  assigned  to  another  grid  point. 

Even  in  this  case,  the  observation  will  always  be  moved  less  thai  one  grid  interval 
(see  Figure  4-1) : 


oX, 
/  3 
#3 


Figure  4-1.  Nearest  Neighbor  Analysis.  "X's"  are  observations  on  the  indicated 
grid  of  points.  is  the  closest  to  Point  #2  but  #2  is  already  assigned  to  a 
closer  observation  X2  so  X^  is  assigied  to  the  next  nearest  point,  # 3. 

(b)  Step  2.  Fill  in  the  values  at  unassigned  points  using  linear  interpolation  of  the 
assigned  values  at  surrounding  points.  The  linear  interpolation  becomes  a  simple 
averaging  at  regularly  spaced  grid  points.  Figures  4-2(a)  and  (b)  are  examples  of 
two  point  and  four  point  averaging  to  fill  in  a  midpoint. 
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Figure  4-2(a) . 

Four  Point  Fill-In 


Figure  4-2(b). 

Two  Point  Fill-In 


When  a  point  has  been  "filled  in,  it  can  then  be  used  as  a  value  to  assist 
in  filling  in  of  other  adjacent  inassijjied  points.  Thus,  more  that  one  pass 
through  the  data  may  be  required  to  fill  in  all  the  points.  For  points  where 
interpolation  is  not  possible,  grid  points  can  be  assigned  a  value  which  is 
the  average  of  grid  point  values  over  the  set  of  grid  points  that  have  already 
been  assigned  values  and  are  located  exactly  one  grid  interval  from  the  grid 
point  of  interest. 

(c)  Step  3.  Smooth  the  grid  point  values  using  a  smoothing  operator  as  given  by 
Equation  4-1-  Suitable  smoothers  can  be  used  at  grid  boundary  and  comer 
points  (e.g. ,  3  or  2  point  smoothers) . 

S  *  (4  x  G,  +  Gj  +  CL.  +  4  +  G4)/*»  4-1 

-K3, 

X 

C.+  -K  +G^ 

-  c 

+C., 

Figure  4-3.  Grid  noint  values  used  in  smoothing  G  point. 

Tht  three  stej  procedure  is  the  same  for  both  surface  and  upper  air  analyses  except  for  a  modification  to 
Step  1  for  upper  air  analyses.  The  value  from  the  upper  ail  observation  is  assigned  to  the  nearest  grid 
point  as  ir.  the  surface  analysis.  However  data  are  then  expanded  to  selected  surrounding  grid  points. 
The  expanded  number  of  grid  points  with  data  values  assigned  as.  a  result  of  wind  information  shall  be 
based  or.  the  rules  contained  in  Western  Region  Technical  Attachment  No.  32-14  and  are  as  follows: 

1.  No  wind  or  missing  wind  -  no  data  expansion. 

2.  Data  are  expanded  plus  and  minus  two  grid  points  in  the  X-  and  Y- direct ions  if: 

(a)  Wind  direction  (with  respect  to  the  grid)  is  300-330'  .  030-060°. 

120-150°,  210-240°,  regardless  of  speed  (Figure  4-4) . 

(b)  Speed  is  less  than  10  m/sec  above  700  nib  and  less  than  5  m/s ec  a" 

700  mb  and  less  thai-  5  m/sec  a'  700>  and  Dele*'  regardless  *nriu 
direction  (Figure  4-4) . 

2.  Data  are  expanded  three  points  in  the  X-direction  and  one  in  the  Y-directian  if 
the  wind  direction  is  240- 30^°  or  060-120°  ano  the  win  sneec  :.s  greater'  than 
10  m/'sec.  above  700  nib  and  greater  than  5  m/sec  at  or  Sclav  700  mb  (Figure  4-5). 

4  Data  are  expanded  three  points  in  the  Y- direct ion  and  env  in  the  X-direction  if 
the  wind  direction  is  330-030  or  150-210'  and  the  wind  speed  is  as  specified  in 
3,  above  (Figure  4-C  . 
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Figure  4-6.  Data  Expansion  from  +  When  Winds  are  North  or  South. 

The.  value  assigned  to  all  the  expanded  points  is  the  same  as  the  central  point,  for  all  fields  except 

height  of  the  constant  pressure  surface.  For  heights,  a  geostrophic  gradient  can  be  calculated  frcni  the 
observed  wind  and  this  neight  value  is  used  to  modify  the  values  that  are  assigned  to  the  surrounding 
ooints.  Tie  basic  equation  for  computing  the  geostrophic  gradient  is 

32/  3n  -  c  x  (f/g)  4-2 

where  3Z/3n  represents  the  difference  in  height  over  distance  n  (perpendicular  to  the  wind  direction; 
i.e.,  height  gradient,  c  is  the  wind  speed,  f  is  the  coriolis  parameter,  €  =  2  a  sin  b  where  w  =  '.29 
l',- 5  anc.  fc  is  latitude,  and  g  is  the  acceleration  of  gravity  g  =  9.806  m/s'" 

To  assign  height  values  to  expanded  points  on  the  grid,  it  is  necessary  to  obtain  u  anc.  v  conpooents  o, 

tne  wind  (with  respect  to  the  grid)  from  wind  direction  and  wind  speed  using  tne  appropriate  equation,;  to: 

the  grid  ix  use  from  Section  2.1.1.  The  value  of  height  Z  at  point  (I+i,  J+j)  where  i  and  j  are  positive 

or  negative  integer  grid  spacing  increments  is  then  computed  from  the  value  at  height  2.  at  point  .J" 

using  tne  following  equation .  p®r~ 


26 


4-3 


Z(I+i,  Jfj)  =  Z(I,J)  +  A  Z  +  A  Z 

x  y 


where 


A  Zx  =  bz  i  dg  v  f/g  4-4a 

A  Zy  =  by  j  de  u  f/g  4 -4b 

where  bx  and  b  are  given  in  Table  4-1  for  each  type  of  grid  and  d  is  the  grid  spacing  at  point  (I, J) . 
The  grid  spacing  dg  is  variable  and  is  a  function  of  the  particular  grid  and  the  location  within  the  grid. 
See  the  discussion  below  Equation  2-11  in  Section  2.2.2  for  the  method  to  calculate  d  . 

If  an  expanded  grid  point  falls  on  cue  or  more  other  expanded  grid  points,  the  average  of  the  two  or  more 
expanded  grid  point  values  should  be  used.  If  an  expanded  grid  point  falls  on  an  original  grid  point  (a 
grid  point  assigned  via  Step  1),  the  original  grid  point  value  is  used. 

The  nearest  neighbor  procedure  can  be  applied  using  a  grid  interval  equivalent  to  that  of  the  SGDB  half¬ 
mesh  polar  stereographic  grid  and  the  SGDB  half-mesh  tropical  grid. 

Two  smoothing  passes  may  be  made  using  the  smoothing  operator,  Equation  4-1. 

Table  4-1 

Values  of  b  and  b 
x  y 

Grid  Type  bx  by 

Northern  Hemisphere  Polar  +1  +1 

Stereographic  Grid 

Southern  Hemisphere  Polar  +1  +1 

Stereographic  Grid 

Mercator  SGDB  Tropical  -1  -1 

4.2.2  Geos trophic  Winds  and  Vorticity.  In  some  cases  it  may  be  necessary  to  derive  the  wind  field 
from  the  height  field  on  a  constant  pressure  surface.  A  conmon  approximation  for  this  purpose  is  the 
geos trophic  approximation,  in  which  it  is  assuned  there  is  exact  balance  between  the  pressure  gradient 
force  and  the  coriolis  force.  The  geostrophic  winds  are  confuted  by  inverting  Equations  4-4a  and  4-4b; 
that  is 


vg  =  A  Zx  g/(f  bx  i  dg)  4 -4a* 

ug  =  A  Zy  g/(f  by  j  dg)  4-4b' 


where  the  subscript  g  assigned  to  u  and  v  indicates  these  are  geostrophic  wind  components. 


The  geostrophic  vorticity  (  t  )  j.s  a  measure  of  the  area-integrated  circulation  of  the  atmosphere.  In  the 

O 

Northern  Hemisphere  c  is  positive  in  lew-pressure  systems  and  troughs  (cyclonic  circulation)  and  is 
S 

negative  in  hig^i-pressure  systems  and  ridges  (anticyclonic  circulation) .  It  is  computed  using  the 
equation: 


where  Av  i  (A  u  )  is  the  difference  of  v  (u  )  between  the  points  I  +  1  (J  +  1)  and  I  -1  (J  -1)  and  c 
§  S  §  §  S 

applies  at  the  point  I,J. 

A. 2. 3  Vertical  Velocity.  The  atmospheric  vertical  velocity  (w)  is  one  of  the  most  important  dynamic 
variables  for  diagnosing  and  predicting  atmospheric  processes.  Unfortunately,  the  vertical  velocity  can¬ 
not  be  measured  directly  on  a  large  scale  with  presently  available  instruments.  Rather,  it  must  be 
inferred  or  computed  from  other  variables  which  are  measured,  such  as  horizontal  winds  or  temperatures. 
There  are  basically  three  methods  for  computing  the  vertical  velocity: 

a.  The  adiabatic  method.  This  method  is  based  on  the  first  law  of  thermodynamics  under  the 
assunption  of  isentropic  flow.  It  requires  computing  time  derivatives  along  with  horizontal  advection. 

As  radiosonde  data  are  normally  available  only  each  12  hours,  there  are  serious  questions  regarding  the 
applications  of  this  method.  Thus,  use  of  the  adiabatic  method  is  discouraged. 

b.  The  omega-equation.  This  method  is  based  on  the  full  set  of  dynamical  equations,  often 
under  the  assunption  of  isentropic  flow.  It  has  the  advantage  of  requiring  observations  only  at  one 
time;  i.e.,  no  time  derivatives  need  to  be  computed.  However,  its  solution  requires  inverting  an  operator 
similar  to  the  three-dimensional  Laplacian.  While  this  can  be  done,  (e.g.,  as  in  NOAA-TM-NWS  WR-138, 
February  1979),  it  is  a  prodigious  task  requiring  the  proper  application  of  boundary  conditions.  Compari¬ 
son  studies  in  the  scientific  literature  find  no  certain  evidence  that  the  omega-equation  produces  irore 
reliable  estimates  of  vertical  velocity  than  other  methods  (e.g.,  Staith  and  Lin,  1978,  Fbnthly  Weather  Re¬ 
view,  Volume  106,  pp.  1687-1694) .  In  view  of  these  faots,  use  of  the  omega  equation  by  local  units  is 
discouraged. 

c.  The  Kinematic  method.  This  method  is  based  on  the  equation  of  continuity;  i.e.,  mass  is 
neither  created  nor  destroyed,  and  is  the  most  straightforward  method.  The  success  or  failure  of  this 
method  to  give  reliable  results  depends  primarily  on  how  well  the  observations  of  horizontal  winds  repre¬ 
sent  actual  atmospheric  conditions.  Very- small-scale  features  of  the  horizontal  winds  must  be  reduced  by 
applying  a  smoothing  process  (as  in  Equation  4-1)  before  computing  vertical  velocities.  The  Kinematic 
method  can  be  applied  to  estimate  vertical  velocities  at  all  heights  once  the  horizontal  wind  fields  have 
been  gridded  and  smoothed  as  described  in  Section  4.2.1.  The  first  step  is  to  estimate  the  vertical 
velocity  at  the  top  of  the  planetary  boundary  layer  (wo>  in  m/s)  according  to: 

wq  -  (K/2f)%  c0  +  (u*(Ax  He)  +  v*(Ay  H0)}/de  4-6 

9 

vdiere  K  =  10m  /s,  f  is  the  coriolis  parameter,  cq  is  the  vorticity  at  the  first  level  above  the 

surface,  u*  and  v*  are  composite  winds  which  have  the  sane  speed  as  the  surface  wind  but  vhich  have  the 

same  direction  as  the  winds  2000  feet  above  the  surface.  A„  H  and  A  H  are  the  slopes  of  the  terrain 

x  e  y  e  r 

along  the  *  and  y  axes  of  the  grid  system  at  the  station.  These  slopes  mist  be  taken  over  spatial  scales 
the  sane  size  as  the  spacing  between  points  of  the  grid  system  (i.e.,  over  distance  dg) .  The  vertical 
velocity  midway  between  layer  k  and  layer  kfl  is  computed  frem 
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Vl'\'AZk{  (AX  u/bx  1  de>  +  <Ay  v/by  j  de>}  (4-7)  j 

where  Axu  is  the  difference  of  u  between  grid  points  on  either  side  of  the  station  along  the  x  axis,  and 
Ayv  is  the  difference  taken  along  the  y  axis.  is  one-half  the  difference  in  height  between  layers 
k+1  and  k-1.  The  process  is  illustrated  in  Figure  4-7; 


layer  k  =2 

the  computed  w  applies  to 
the  altitude  midway  between 
layer  k  =  1  and  k  =  2. 

Horizontal  - >  .  < - layer  k  =  1 

convergence 

layer  k  =  0 

Figure  4-7.  Illustration  of  calculation  of  vertical  velocity. 

The  vertical  velocity  at  successively  higher  layers  is  computed  by  resolving  Equation  4-7.  A  note  of 
caution  must  be  given.  The  vertical  velocities  computed  by  the  above  algorithm  are  useful  for  qualitative 
analysis  because  the  algorithm  given  here  is  cast  in  height  rather  than  pressure  coordinates.  That  is,  in 
the  lower  and  mid- troposphere  the  sign  of  the  vertical  velocity  should  be  reliable,  except  for  very  small 
values,  and  the  magnitude  of  the  motion  should  permit  description  as  "small"  or  "large."  Finally,  as  with 
any  other  analysis  scheme,  meteorologists  will  appreciate  the  uses  and  shortcomings  of  this  variable  only 
when  applied  in  the  context  of  a  complete  meteorolgical  analysis. 

4.3  Cross  Section  Objective  Analysis.  Vertical  cross  section  analysis  is  a  tool  used  by  meteorolo¬ 
gists  for  depicting  the  spatial  and  temporal  distribution  of  temperature,  moisture,  and  wind  velocity. 
Radiosonde  data  obtained  at  upper  air  observation  stations  are  input.  The  cross  section  is  very  useful  in 
that  it  shows  subsynoptic  scale  gradients  from  the  detailed  vertical  information  at  the  upper  air  stations 
which  are  synoptical ly  spaced. 

4.3.1  Hermite  Polynomial  Interpolation  for  Objective  Cross  Section  Analysis.  One  of  the  most  useful 
types  of  cross  sections  is  one  that  depicts  the  potential  temperature,  0.  It  is  called  an  isentropic 
section  and  it  has  the  characteristic  of  showing  meteorological  regions  of  interest  such  as  frontal  zones, 
inversions  and  the  tropopause  in  significant  detail.  Near  these  regions,  there  are  strong  gradients  of  9. 
The  one  analysis  method  that  is  most  widely  used  because  it  is  both  relatively  simple  and  accurate  is 
Hermitian  interpolation  described  by  Shapiro  and  Hastings  (1973)  in  Journal  of  Applied  Meteorology, 

Volume  12,  pp.  753-762.  To  graph  the  spatial  distribution  of  pressure  along  an  isentropic  surface,  h(x) , 
Hermite  interpolation  polynomials  can  be  used  to  define  the  isentropic  interpolation  function,  k(x).  The 
Hermite  interpolation  polynomials  are  cubic  polynomials  pieced  together  at  selected  points  where  both 
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function  values  and  first  derivatives  are  known.  8  and  0g  are  the  two  parameters  to  be  analyzed  using 
Hermitian  polynomials.  Values  of  0  and  0g  are  determined  from  Equations  3-1  and  3-6  ,  respectively.  All 
equations  used  in  the  Hermite  interpolations  and  the  explanations  for  boundary  conditions  and  differencing 
schemes  are  contained  in  the  Shapiro  and  Hastings  paper. 

4.3.2  Objective  Cross  Section  Analysis  Using  Linear  Interpolation.  Many  fields  in  the  vertical  can, 
at  times,  be  irregular  in  nature.  Many  meteorologists  have  found  that  these  fields  can  De  displayed 
well  on  a  cross  section  using  linear  interpolation.  Interpolation  is  in  two  dimensions,  X  and  log  P. 
Interpolation  in  the  vertical  may  use  standard  and  significant  level  data  to  determine  parameter  values  at 
50  mb  increments  over  a  suitable  range;  e.g. ,  1050  to  100  mb.  The  interpolation  equation  for  quantity  Q  at 
pressure  P  =  (1050  mb  -  n  x  50  mb)  where  n  is  a  integer  such  that  0  _<  n  <  19  is  as  follows : 

Q  =  Qi  +  (Qi+1  -  Qi)  (In  P  -  In  Pi)/(ln  P.+r  In  P.)  4-8 

where  P^  and  Pi+1  are  reporting  levels  adjacent  to  level  P  such  that  >  P  >  P^  Interpolation  in  the 
horizontal  uses  data  from  reporting  stations  or  data  points  computed  from  gridded  data  to  determine  param¬ 
eter  values  at  regularly  spaced  grid  intervals  (half -mesh  grid  spacing)  along  the  horizontal  axis  of  the 
cross  section.  The  interpolation  equation  for  quantity  Q  at  distance  Z  =  n  x  d  along  the  cross  section 
axis  from  the  origin,  where  n  >  0  is  an  integer  and  d  is  the  half -mesh  grid  spacing,  is  as  follows: 

Q  =  Qi  +  (Qi+1  -  V  (x  -  si)/(Xi+1  -  X.)  4-9 

where  X^  and  X^+^  are  distances  to  reporting  stations  or  data  points  computed  from  gridded  data  adjacent  to 
distance  X  such  that  X^  <  X  <  The  parameters  for  which  bi-linear  interpolations  can  be  performed  for 

contouring  on  cross  sections  are  relative  humidity,  dewpoint  depression,  wind  speed,  vertical  wind  shear, 
temperature,  D- values  and  mixing  ratio.  Once  data  are  interpolated  to  (X,P)  grid  points,  a  graphic  repre¬ 
sentation  of  the  parameter  field  can  be  accomplished  using  contours  that  are  unbroken  piecewise  continuous 
between  adjacent  grid  boxes.  One  especially  useful  relation  for  the  analysis  of  vertical  profiles  of  wind 
on  cross  section  representations  is  the  thermal  wind  equation.  The  corrponents  of  this  equation  show  the 
relation  of  horizontal  temperature  gradient  to  vertical  wind  shear.  In  generalized  X,  Y  coordinates,  the 
corrponents  are  given  by: 

3T  =  |T  3V  4- 10a 

3X  g  3Z 

3T  _  _  fT  3U  4-10b 

3Y  g  3H 

where  T  is  usually  estimated  by  the  temperature  midway  on  the  vertical  layer  across  which  the  wind  shear 
is  corrputed. 

5.  mP  PROJECTIONS 

Because  reporting  station  locations  are  given  in  latitude  and  longitude,  plotting  of  parameters  at  the 
geographic  location  of  a  station  on  a  map  requires  a  transformation  to  the  coordinate  systan  used  for  the 
map.  A  required  first  step  in  this  process  is  to  transform  station  locations  from  latitude  and  longitude 
(0,  X)  to  coordinates  (I,  J)  on  the  appropriate  AFX5^3  whole-mesh  Satellite  Global  Data  Base  (SGDB)  grid: 
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the  Northern  Hemispheric  polar  stereographic  grid,  the  Southern  Hemispheric  stereographic  grid,  or  the 
tropical  (Mercator)  grid.  Precise  location  of  stations  on  maps  is  required  and  is  acconplislied  through  the 
use  of  real-valued  SGDB  grid  coordinates.  The  transformation  equations  for  both  the  Northern  and  Southern 
Hemispheric  SGDB  whole-mesh  polar  stereographic  grids  and  for  the  SGDB  whole-mesh  tropical  grid  are  given 
in  AFGMC/TN-79/003.  Units  requiring  access  to  this  TN  should  submit  requests  through  their  parent  unit. 

6.  MILITARY  GRID  REFERENCE  SYSTEM  AND  LATITUEE/LDNGITUDE  COORDINATE  TRANSFORMATIONS 

6.1  General .  The  Military  Grid  Reference  System  (MGRS)  is  a  method  for  identifying  geographical 
points  on  earth.  This  section  describes  the  algorithms  for  converting  latitude/ longitude  to  MGRS  points. 
The  MGRS  uses  the  Universal  Transverse  Mercator  (UIM)  projection.  The  range  of  latitude  it  covers  is  from 
80°  south  to  80°  north.  The  globe  is  divided  into  large  geographical  areas  extending  6°  in  longitude  and 
8°  in  latitude,  each  of  which  is  given  a  unique  identification  called  the  grid  zone  desigiation.  The  6°  by 
8°  areas  are  further  subdivided  into  100,000  meter  squares,  with  each  square  given  a  unique  two-letter 
identification.  Numerical  references  within  these  100,000  meter  squares  are  given  to  100m  resolution  in 
terms  of  "Easting"  and  "Northing"  grid  coordinates  of  the  point. 

6.2  Method  for  Calculating  the  Latitude /longitude  from  MGRS .  Necessary  inputs  to  the  computations 

are: 


.  Applicable  UTM  grid  zone, 

.  Location  (Easting  and  Northing)  of  the  origin  (lower  left  comer)  of  the  pertinent  100,000 
meter  grid  square, 

.  Position  (Easting  and  Northing)  of  the  point  within  the  grid  square  for  the  100m  resolution. 
These  inputs,  that  define  locations  in  the  MGRS,  are  given  in  the  following  format: 

ZZ  LEN  eee  nnn 

where  ZZ  is  the  zone  number, 

L  is  the  8°  latitude  band  in  which  the  point  is  located, 

E  defines  the  100,000m  Easting  division, 

N  defines  the  100,000m  Northing  division, 

eee  further  defines  the  Easting  location  to  the  nearest  100m, 

nnn  further  defines  the  Northing  location  to  the  nearest  100m. 

A  listing  of  variables,  constants,  and  notations  used  in  the  conversion  algorithms  for  the  MGRS: 

Cl  =  ZZ  =  Zone  number 

C2  =  L  =  Eight  degree  latitude  band  (an  alpha  character) 

C3  =  E  =  Easting  division  (an  alpha  character) 
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C4  -  N  »  Northing  division  (an  alpha  character) 

C5  -  eee  -  Easting  location 
C6  *  rmn  -  Northing  location 
LAT  »  Latitude  (in  decimal  degrees) 

LON  =  Longitude  (in  decimal  degrees) 

HO  =  6,378,338 
Rl  -  0.00672267 
R2  =  0.9996 
R3  =  0.048481368 

R4  =  3600 
R5  =  500,000 
R6  =  6367645.45 
R7  *  16106.99 
R8  =  16.976 

R19  =  0.017453295 
R25  =  10,000,000 
R26  =  100,000 
R27  -  100 

**  is  notation  for  "raised  to  the  power  of." 
x  is  notation  for  "multiplied  by." 

/  is  notation  for  "divided  by." 

ABS  is  notation  for  "take  the  absolute  value  of." 

9QRT  is  notation  for  "take  the  square  root  of." 

1ST  is  notation  for  "take  the  integer  of'  (discarding  the  fractional  part  of  a  nunber) . 


Latitude  is  determined  from  Equation  6-1  : 


LAT  =  02'/ 10000  x  INT  ((RIO  -  R12  x  R9  x  R9  -  R13  x  (R9  **  4))/R4)  x  10000  +  .5)  6-1 

R4  is  a  constant 

Each  of  the  terms  in  Equation  6-1  is  determined  from  tables  and  other  equations  according  to  the  follow¬ 
ing  steps: 

Step  1:  To  determine  C2' ,  find  CX  from  Table  6-1.  NOTE:  The  range  of  C4'  for  use  in  Step  3  below. 

Table  6-1  -  CX 


C2 

CX 

C4'  Range 

C2 

CX 

C4'  Raige 

C 

-10 

11-20 

N 

0 

00-08 

D 

-9 

20-28 

P 

1 

08-17 

E 

-8 

28-37 

Q 

2 

17-26 

F 

-7 

37-46 

R 

3 

26-35 

G 

-6 

46-55 

S 

4 

35-44 

H 

-5 

55-64 

T 

5 

44-53 

J 

-4 

64-73 

U 

6 

53-62 

K 

-3 

73-82 

V 

7 

62-71 

L 

-2 

82-91 

W 

8 

71-79 

M 

-1 

91-99 

X 

9 

79-88 

greater 

than  or 

equal  to  zero,  then  C2' 

=  1. 

If  CX  is  less 

than  zero,  then  C2 

Step  2:  Determine  C3'  from  Table  6-2. 


Step  4:  Confute  the  RIO  term  in  Equation  6-1  from  Equation  6-2  : 


r" 


RIO  -  ((R17/R2)  +  R7  x  sin  (2  x  R17)(R6  x  R2  x  R19)) 

-  R8  x  sin  (4  x  R17/(R6  x  R2  x  R19))  x  (1/(R6  x  R19)) 


R2 ,  R7,  R6,  R19,  and  R8  are  constants. 

Step  5:  Compute  the  R17  term  in  Equation  6-2  from  Equation  6- 3a  : 

If  C2'  is  less  than  zero,  then: 

R17  -  R25  -  C6  x  R27  -  C4'  x  R26  6-3a 

R25,  R27,  and  R26  are  constants,  C6  is  the  Northing  location  and  C4‘  was  determined  in  Step  3. 
If  C2'  is  greater  than  or  equal  to  zero,  then: 

R17  =  C6  x  R27+  C4'  x  R26  6-3b 

Step  6:  Compute  the  R12  term  in  Equation  6-1  from  Equation  6-4  : 

R12  =  R16  X  (10  **  16)/  (2XR3XR2XR2X  Rll  X  Rll)  6-4 

where:  R2  and  R3  are  constants,  and  R16  and  Rll  are  determined  from: 

R16  -  TAN  (RIO)  6-5 

Rll  =  RO/SQRT  (1-R1  x  sin  (R17/(R6  x  R2  x  R19))  x  sin 

(R17/(R6  x  R2  x  R19)))  6-6 

RO,  Rl,  R6,  R2,  and  R19  are  constants  and  R17  was  previously  defined  (Equations  6- 3a  or  6- 3b 
above) 

Step  7:  Compute  R9  in  Equation  6-1  from  Equation  6-7  : 

R9  =  (R18  -  R5)  x  (0.000001)  6-7 

tdiere  R18  is  given  by: 

R18  =  C5  x  R27  +  C3'  x  R26  6-8 

C5  is  the  Easting  location,  R27  and  R26  are  constants  and  C3'  was  determined  in  Step  2. 

Step  8:  Conpute  R13  in  Equation  6-1  from  Equation  6-9  • 

R13  -  (R16  X  R16  X  3  +  5)  X  R16  X  (10  **  28)/ 

(24  X  (R2  **  4)  X  R3  X  (Rll  **  4))  6-9 
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R2  and  R3  are  constants,  R16  and  Rll  are  defined  by  Equations  6-5  and  6-6  ,  respectively. 

If  latitude  in  Equation  6-1  is  negative,  it  is  south  latitude. 

Longitude  is  determined  from  Equation  6-10  •. 

LON  -  (1/10000)  x  LOT  (((R9  x  (R14  -  R15  x  R9  x  R9)/R4) 

+  R22)  x  10000  +  0.5)  6-10 

If  ION  is  negative,  it  is  west  longitude. 

Step  1:  R9  is  given  by  Equation  6-7. 

Step  2:  R14  is  determined  from  Equation  6-11; 

R14  =  (10  **  10) /(R2  x  R3  x  Rll  x  cos  (R10))  6-11 

R2  and  R3  are  constants,  Rll  is  given  by  Equation  6-6  and  R10  is  given  by  Equation  6-2. 

Step  3:  R15  in  Equation  6-10  is  determined  from: 

R15  =  (R16  x  2  +  1)  x  (10  **  22)/ (6  x  R3  x  R2  x  R2  x 

R2  x  Rll  x  Rll  x  Rll  x  R23)  6-12 

R16  is  given  by  Equation  6-5  ,  R2  and  R3  are  constants  and  Rll  is  defined  by  Equation  6-6.  R23  i 

determined  from-. 

R23  =  cos  (R10)  6-13 

Step  4:  R22  in  Equation  6-10  is  determined  from: 

R22  =  Cl  x  6  -  183  6-14 

6.3  Method  for  Calculating  the  MSRS  Location  from  Lat itu de/longitude .  To  determine  the  MGRS  point, 


calculate  Cl.  C2,  C3,  C4.  C5,  and  C6. 

Cl  =  R21  6-15 

where: 

R21  =  INT  (LON/6  +  31)  6-16 

The  following  R's  will  be  used  in  the  determinations  of  C2  through  C6.  (These  values  of  R  are  new  sets  of 
R's  from  those  in  Section  6.2,  although  the  constants  given  in  that  section  remain  the  same.) 

R9  -  -((R21  x  6  -  183  -  LON)  x  0.36  6-17 


R10  -  LAT,  inless  LAT  is  less  than  zero,  then  R10  -  -LAT 


Rll  «  RD/SqKT  (1-R1  x  sin  (RIO)  x  sin  (RIO)) 

R12  =  R2  x  (RIO  x  R19'  x  R6  -  R7  x  sin  (2  x  RIO)  + 

R8  x  sin  (4  x  RIO)) 

R13  =  (R2  x  R3  x  R3  x  Rll  x  sin  (2  x  R10))/4 

R14  =  R2  x  (R3  **  4)  x  (1/24)  x  Rll  x  sin  (RIO)  x  cos  (RIO) 
x  cos  (RIO)  x  cos  (RIO)  x  (5  -  tan  (RIO)  x  tan  (RIO)) 

R15  =  R2  x  R3  x  Rll  x  cos  (RIO) 

R16  -  R2  x  (R3  **  3)  x  (1/6)  x  Rll  x  cos  (RIO) 
x  (2  x  cos  (RIO)  x  cos  (RIO)  -  1) 

R17  =  (R14  x  (R9  **  4)  +  R13  x  R9  x  R9  +  R12) 
x  LAT/ABS  (LAT) 

R29  =  ETC  (((R9  **  3)  x  R16  +  R15  x  R9  +  R5)/R26) 

R18  »  (R9  **  3)  x  R16  +  R15  X  R9  +  R5  -  R29  x  R26 

If  R17  is  less  than  zero,  then: 

R  =  R17  +  R25 

C2'  =  INT  (-1  -RIO/ 8) 

If  R17  is  greater  than  or  equal  to  zero,  then: 

R  =  R17 

C2'  *  INT  (R10/8) 

To  determine  C3' : 

If  (R21  +  2)/3  =  INT  ((R21  +  2)/3)  then  C3'  -  R29  +  10 

If  (R21  +  1  )T-  =  INT  ((R21  +  l)/3)  then  C3’  =  R29  +  20 

Otherwise,  C3'  =  R29  +  30 
C4'  is  given  by: 

C4'  =  INT  (R/R26) 

C5  and  C6  are  determined  from  Equations  6-36  and  6-37  ,  respectively: 


C5  *  INT  (R18/R27) 


Co  =  INT  ((R  -  04'  x  R26)/R27) 
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At  this  point,  Cl,  C2’,  C3’ ,  C4‘ ,  C5,  and  C6  are  taown. 


C2  is  determined  from  Table  6-4: 


Table  6-4  -  C2  and  C2’ 


C2 


C2’ 


C 

D 

E 


G 

H 


-10 

-9 

-8 

-7 

-6 

-5 


-3 


M 


-1 


Determine  C3  from  Table  6-5  where  C3’  =  A//B 
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7.  SKEW-T,  LOG-P  BACKGROUND  LINES  AND  TEMPERATURE  AND  DEWPOINT  TEMPERATURE  PROFILES 


7.1  General .  The  Skew-T,  Log-P  diagram  is  discussed  in  AWS/TN- 79/006.  It  is  used  for  the  graphical 
depiction  of  the  vertical  temperature  and  moisture  (dewpoint  temperature)  structure  of  the  atmosphere  over 
a  reporting  station.  It  is  also  used  for  plot  and  display  of  winds  at  standard  constant  pressure  surfaces 
and  heights  of  these  surfaces  above  mean  sea  level.  Derived  parameters  as  in  Chapter  3  can  also  be  display¬ 
ed  on  the  Skew-T,  Log-P  diagram.  To  produce  the  graphical  depiction  on  a  computer  display  terminal  it  is 
necessary  to  develop  a  version  of  the  Skew-T,  Log-P  background  diagram  comprised  of  four  sets  of  lines 
upon  which  the  vertical  profiles,  plots,  and  parameter  values  will  be  displayed.  The  equations  for  the 
background  lines  and  for  generating  profiles  of  temperature  and  dewpoint  temperature  are  contained  in  the 
following  subsections.  The  material  in  these  subsections  has  either  been  taken  directly  (quotes)  or  ab¬ 
stracted  from  the  document  "Algorithms  for  Generation  of  a  Skew-T,  Log-P  Diagram  and  Computing  Selected 
Meteorological  Quantities,"  by  G.  S.  Stipanuk  (1973)  and  published  by  the  Army  Electronics  Conmand  (E00M- 
5515) . 

7.2  Background  Lines  and  Vertical  Profiles  for  Temperature  and  Dewpoint  Temperature.  Four  sets  of 
lines  comprise  the  background  Skew-T,  Log-P  diagram  —  pressure,  temperature,  dry  adiabat  and  saturation 
adiabat  lines.  The  first  two  sets  of  curves,  temperature  and  pressure,  are  used  to  locate  points  on  the 
diagram.  An  arbitrary  coordinate  systen  is  used  to  measure  distances.  The  origin  corresponds  to  the  point 
at  a  temperature  of  0  degrees  Celsius  (°C)  and  a  pressure  of  1000  mb.  The  X-direction  is  parallel  to  the 
pressure  lines  (horizontal)  with  positive  X  to  the  right.  The  Y-direction  is  perpendicular  to  the  X- 
direction  with  positive  Y  towards  lower  pressures  (up) .  A  point  cn  the  diagram  specified  by  temperature 
and  pressure  is  transformed  to  X,  Y  coordinates  by  Equations  7-1  and  7-2  : 

X  =  . 1408T  -  10.53975  log  P  +  31.61923  7-1 

Y  =  -11.5  log  P  +  34.5  7-2 

The  components  in  the  X,  Y  coordinate  system  are  given  in  inches.  For  most  applications,  the  coordinate 
system  for  pressure  ranges  from  1050  to  100  mb  and  for  tenperature  it  ranges  frcm  -40°C  to  440°C  along  the 
1000  mb  pressure  line.  Because  the  temperature  lines  are  skewed,  this  30°  Celsius  range  will  slide  to 
lower  temperatures  for  decreasing  pressure  values.  It  will  be  necessary  to  scale  the  dimensions  given  by 
Equations  7-1  and  7-2  such  that  the  entire  background  Skew-T,  Log-P  diagram  (1050  to  100  mb  and  440°  to 
-40°C  at  1000  mb)  can  be  displayed  within  the  dimensions  of  the  graphic  display  device.  The  equations  for 
dry  adiabats  and  for  saturation  adiabats  are  given  in  Table  7-1.  The  tenperature  T  at  an  arbitrary 
pressure  on  a  saturation  adiabat  is  determined  by  the  bisection  method.  The  bisection  method  is  a  numeri¬ 
cal  technique  which  decreases  the  difference  between  the  upper  and  lower  estimates  by  a  factor  of  %  per 
iteration.  The  temperature  is  assured  to  lie  in  the  range  -80°C  to  40°C.  An  initial  guess  of  -20°C  is 
made  and  the  correction,  T*,  computed.  The  correction  term  decreases  by  a  factor  of  \  after  each 
correction.  Terminating  after  13  corrections  usually  gives  satisfactory  results.  The  algorithm  for  com¬ 
puting  the  temperature  on  a  saturation  adiabat  is  based  on  Equation  7-3  : 

0  =  (OJ  x  eyp  (-Lr  /C_T)  7-3 

6  S 

In  addition  to  the  algorithms  which  generate  the  curves  for  each  family,  it  is  necessary  to  have 
algorithms  which  determine  which  curve  in  a  family  passes  through  an  arbitrary  point  (T,  P).  Algorithms  to 
accomplish  this  are  given  in  Table  7-2. 


SKEW-T,  LOG-P  Algorithms 


Fanily 

Dry 

Adiabat,  (T^ 


Saturation 
Adiabat ,  (T<^) 


Parameter 

Algorithm 

0,  the  potential 

’W9,p)  =  e(p/iooo)0-2854 

tenperature 

T  is  in  °Kelvin.  °K=°C  + 

273.16 

0g,  the  Temperature 

TSA  (6s’  P>  “  T1  +  ^  Tt 

at  1000  nfc 

i=l 

T*  =  120/21  SIGNta  exp[brs  (Tt>P)/Ti] 

-  Tt  (1000/P)0-2854 


T.  =  T .  i  +  T*  7-7 

1  1‘1  i-1 

a  =  0g  Which  is  given  by 
Equation  7 -9 

b  =  -7551.8986 

r  (T.  P)  =  0.622e  /(P-e  ) 
s  s  s 

eg  is  given  by  Equation  C-2(b) 
or  C-2(c) 


The  SIGN  function  is  -1  or  +1  corresponding  to  the  algebraic  sign  of  the  argument. 


Table  7-2 


Determining  a  Curve  Through  a  Given  Point 

Curve  Set  Parameter  for  Curve  Passing 

Through  (T,  P) 


Dry  Adiabat 

6  =  T  (1000/P)0'2854 

7-8 

Saturation  Adiabat 

0s  =  T  (1000/P)0'2854/ 

7-9 

exp  [br.(T,  P)/T] 
b  =  -2651.8986 


8.  INSOLATION  AND  ASTRONOMICAL  ALGORITHMS 


a.  Hour  Angles.  The  following  formulas  are  useful  if  astronomical  data,  such  as  that  given 
Navigational  and  Stellar  Tables  are  applied  to  navigational  purposes: 

GHA  =  15(GAST  -  RA) 


LHA  =  15CLAST  -  RA)  =  (HA  -  X 


GHA  Aries  =  15  GAST 
SHA  =  360°  -  15  RA 
GHA  =  (HA  Aries  +  SHA 


where 

GHA  is  the  Greenwich  hour  angle  in  degrees; 

LHA  is  the  local  hour  angle  in  degrees; 

(HA  Aries  is  the  Greenwich  hour  angle  of  the  First  Point  of  Aries  (the  origin  of  right 
ascension)  in  degrees; 

SHA  is  the  sidereal  hour  angle  in  degrees; 

RA  is  the  apparent  right  ascension  (referred  to  the  true  equator  and  equinox  of  date) 

in  hours; 


X  is  the  local  longitude  in  degrees  (west  is  positive;  east  is  negative); 


GAST  is  the  Greenwich  apparent  sidereal  time  in  hours; 


LAST  is  the  local  apparent  sidereal  time  in  hours; 

W,A  are  constants  from  a  tabulated  list. 

The  GHA  and  declination  of  the  sun  and  mocn  are  computed  from  a  power  series  of  the  form 
F(X)  =  aQ  +  ax  X  +  a2  X2  +  a3  X3  +  a4  X4  +  a$  £ 

where  X  is  a  time-like  variable  such  that  -1  <  X  <1  and  the  coefficients  a  to  ac  are  tabulated  constants. 
The  coefficients  change  from  month- to-iaonth  within  each  year,  and  also  change  from  year-to-year.  Uiits 
requiring  the  current  year's  tabulated  coefficients  should  submit  requests  through  channels  to  their  parent 
unit.  To  evaluate  the  power  series  for  one  of  the  navigational  variables,  it  is  first  necessary  to  evalu¬ 
ate  X: 


(1)  Ccmpute  t,  where  t  =  N  +  CM?/ 24  N  is  the  day  of  the  year  at  Greenwich  and  CM  is 

Greenwich  Mean  Time  in  hours .  N  is  tabulated  on  many  calendars ,  or  can  be  confuted 
using  Equation  8-3. 

t- W 

(2)  Compute  X,  where  X  =  (-j~j-l  .  W  and  A  are  constants  tabulated  along  with  the 

coefficients  aQ  to  a^  in  the  Almanac  for  Computers . 

As  an  exanple,  find  the  sun's  GHA  at  19*102m30S  CM  on  23  May  1983:  From  Equation  8-3  we  find  that  N  =  143 

19  0417 

on  23  May.  Thus,  t  =  143  +  — ^ —  =  143.7934.  Next,  a  page  from  the  Almanac  for  Computers  is  shown  in 
Figure  8-1.  Frcm  this  page,  the  appropriate  values  of  W  and  A  are  121  and  16,  respectively.  Thus, 

„  _  143.7923  -  121 
x  ^  -  l  =  0.4145875 

Also,  frcm  Figure  8-1,  the  values  of  aQ  to  a^  for  OTA  for  23  May  1983  are  found  and  entered  in  the  power 
series.  That  is, 


OTA  =  5940.9204  +  X  (5759.9128)  +  X2  (-0.2935)+  X3  (0.0195)  +  X4  (0.0046)  +  X5 
(-0.0041)  =  8386.45605 

This  value  lies  outside  the  range  0°  <  GHA  <  360°,  so  an  integer  multiple  of  360°  is  subtracted  from  OTA. 
In  this  case  the  integer  is  23,  so 

OTA  =  8386.45605  -23  (360)  =  106°. 45605. 


Days  121  thru  152  JD  2445455.5  to  2445487.5  Dates  May  1  thru  June  1 

A  =  16.0  B  =  -8.56250000  W  =  121 


Aries 

Sun 

Sun 

Sun 

GHA 

GHA 

Dec 

S  D 

Term 

O 

O 

O 

O 

0 

5994.1267 

5940.9204 

19.1483 

0.2642 

1 

5775.7719 

5759.9128 

3.6747 

-0.0033 

2 

-0.0002 

-0.2935 

-0.6950 

0.0 

3 

-0.0028 

0.0195 

-0.0522 

0.0061 

4 

0.0001 

0.0046 

0.0064 

0.0 

5 

0.0013 

-0.0041 

0.0010 

-0.0037 

Sums 

11769.8970 

11700.5597 

22.0832 

0.2633 

Days  152  thru  183 

JD  2445586.5  to  2445518.5 

Dates  June  1  thru  July  2 

A 

=  16.0 

B  =  -10.50000000 

W  =  152 

Aries 

Sun 

Sun 

Sun 

GHA 

GHA 

Dec 

S  D 

Term 

O 

O 

O 

O 

0 

6024.6817 

5939.8423 

23.3565 

0.2633 

1 

5775.7718 

5759.1351 

0.5474 

0.0 

2 

0.0009 

-0.0345 

-0.8805 

0.0 

3 

-0.0035 

0.0765 

-0.0072 

0.0 

4 

-0.0010 

-0.0054 

0.0091 

0.0 

5 

0.0023 

-0.0074 

-0.0005 

0.0 

Sums 

11800.4522 

11699.0066 

23.0248 

0.2633 

Days  182  thru  213 

JD  2445516.5  to  2445548.5 

Dates  July  1  thru  Aug  1 

A 

=  16.0 

B  =  -12.37500000 

W  =  182 

Aries 

Sun 

Sun 

Sun 

GHA 

GHA 

Dec 

S  D 

Term 

O 

O 

O 

O 

0 

6054.2507 

5938.4969 

21.3321 

0.2633 

1 

5775.7703 

5759.6299 

-2.6376 

0.0 

2 

0.0027 

0.2866 

-0.7707 

0.0 

3 

-0.0004 

0.0377 

0.0375 

0.0 

4 

-0.0019 

-0.0158 

0.0033 

0.0 

5 

0.0002 

-0.0027 

0.0036 

0.0 

Sums 

11830.0216 

11698.4326 

17.9682 

0.2633 

Days  213  thru  244 

JD  2445547.5  to  2445579.5 

Dates  Aug  1  thru  Sep  1 

A 

=  16.0 

B  =  -14.31250000 

W  =  213 

Aries 

Sun 

Sun 

Sun 

GHA 

GHA 

Dec 

S  D 

Term 

O 

O 

O 

O 

0 

6084.8056 

5938.9411 

13.6746 

0.2633 

1 

5775.7708 

5760.8264 

-5.0615 

-0.0007 

2 

0.0026 

0.2837 

-0.4584 

0.0004 

3 

-0.0014 

-0.0319 

0.0506 

0.0031 

4 

-0.0017 

-0.0100 

-0.0033 

0.0006 

5 

0.0008 

0.0013 

0.0062 

-0.0015 

Sums 

11860.5767 

11700.0106 

8.2082 

0.2652 

Figure  8-1.  Power  Series  Approximation  of  Nautical  Almanac  Data  for  Year  1983. 

When  using  the  above  formulas,  it  may  be  necessary  to  add  or  subtract  360°  to  reduce  the  resulting  hour 
angles  to  the  range  0°  -360 P.  Often  the  local  hour  angle  values  are  reduced  to  the  range  -18CP  to  +180°, 
in  vrtiich  case  they  are  called  meridian  angles.  In  all  cases,  positive  hour  angle  values  are  measured  west 
ward  from  the  meridian. 


b.  Altitude  and  Azimuth.  Notation: 


a  =  altitude  of  body  above  (if  sin  a  >  0)  or  below  (if  sin  a  <  0)  the  horizon; 

A  =  azimuth  of  body  measured  eastward  from  north  over  the  range  CP  <  A  <  360°; 

0  =  latitude  of  observer  (north  is  positive;  south  is  negative); 

6  =  declination  of  body  (north  is  positive;  south  is  negative) ; 

LHA  =  local  hour  angle  of  body; 

z  =  zenith  distance  of  body  (z  =  90°  -  a) . 

The  following  formulas  can  be  used  to  compute  the  altitude  (a)  and  azinuth  (A)  of  a  celestial  body: 

sin  a  =  cos  z  =  sin  0  sin  6  +  cos  0  cos  5  cos  LHA  8-1 

x  =  tan  A  =  sin  LHA/ (cos  LHA  sin  0  -  tan  &  cos  0)  8-2 

Since  computers  and  calculators  normally  give  the  arctangent  in  the  range  -90°  to  +90®,  the  correct 
quadrant  for  A  can  be  selected  according  to  the  following  rules: 

If  0°  <  LHA  <  180°, 

A  =  180°  +  arctan  x,  if  x  is  positive, 

A  =  360°  +  arctan  x,  if  x  is  negative. 

If  180°  <  LHA  <  360°, 

A  =  arctan  x,  if  x  is  positive, 

A  =  130°  +  arctan  x,  if  x  is  negative. 

In  standard  navigational  notation  altitude  and  azimuth  are  denoted  He  and  Zn,  respectively.  Equations 
8-1  and  8-2  are  the  basic  formulas  used  in  preparing  sight  reduction  tables;  they  do  not  include  the 
effect  of  refraction.  EXAMTLE:  Compute  the  altitude  and  azimuth  of  the  Sun  at  19^102m30s  Off  on  23  May 
1983  at  College  Park,  Maryland: 

Latitude:  0  =+39.0®  sin  0  =  +0.6293  cos  0  =  0.77715 

longitude:  X  =  +77.0® 

Using  the  power  series  discussed  earlier,  the  Sun's  GHA  and  6  are  found  to  be 
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hence 


LHA  -  106.457°  -  77.0°  -  29.457° 
cos  LHA  -  +0.87073 


GHA  -  106.457° 

sin  LHA  -  +0.49177 

6  *  +20.579°  sin  6  =  +0.35150  cos  6  -  +0.93619  tan  6  =  +0.37546 

sin  a  -  cos  z  -  (0.62932) (0.35150)  +  (0.77715)(0.93619)(0.87073)  =  +0.85471 

a  -  58.7° 

x  =  tan  A  -  0.49177/((0.87073)(0. 62932)  -  (0.37546) (0.77715) 

-+  1.91963  arctan  x  -  62.5° 

Since  lift  is  greater  thai  0°  and  less  than  180P,  and  since  x  is  positive,  A  -  180°  +  62.5°  =  242.5° 

c.  Sunrise,  Suiset,  and  Twilight.  For  locations  between  latitudes  65°  North  and  65°  South,  the 
following  algorithm  provides  times  of  sunrise,  sixiset,  and  twiligjnt  to  an  accuracy  of  +2m,  for  any  date  in 
the  latter  half  of  the  twentieth  century.  Because  the  phenomena  depend  on  local  meteorological  conditions, 
attempts  to  attain  higher  accuracy  are  seldom  justified.  Although  the  algorithm  can  be  used  at  higher 
latitudes,  its  accuracy  deteriorates  near  dates  on  which  the  Son  remains  above  or  below  the  horizon  for 
more  than  twenty- four  hours.  Notation: 

0  »  latitude  of  observer  (north  is  positive ;  south  is  negative) ; 

A  -  longitude  of  observer  (west  is  positive;  east  is  negative) ; 

M  -  Son's  mean  anomaly; 

L  -  Sun's  true  geocentric  longitude; 

RA  =  Sun's  right  ascension; 

6  =  Sun's  declination; 

H  -  Sun's  local  hour  angle; 

z  =  Sun's  zenith  distance  at  rise,  set,  or  twilight*; 

t  -  approximate  time  of  phenomenon  in  days  since  31  Dec. .c/1  UT; 

T  -  local  mean  time  of  phenomenon; 

UT  -  universal  time  of  phenomenon. 


*Ihe  proper  value  of  z  should  be  chosen  from  the  following: 


z 

cos  z 

Srnrise  and  Sunset 

90.833° 

-0.01454 

Civil  Twilight 

96° 

-0.10453 

Nautical  Twilight 

102° 

-0.20791 

Astronomical  IWi light 

108° 

-0.30902 

Formulas: 

(1)  M  =  0.9856°t  -  3.289° 

(2)  L  =  M  +  1.916°  sin  M  +  0.020°  sin  2M  +  282.634° 

(3)  tan  RA  =  0.91746  tan  L 

(4)  sin  6  =  0.39782  sin  L 

(5)  x  =  cos  H  =  (cos  z  -  sin  6  sin  0)/(cos 6  cos  0) 

(6)  T  =  H  +  RA  -  0.06571ht  -  6.622h 

(7)  UT  =  T  +  A 

Procedure: 

(a)  Step  1.  With  an  initial  value  of  t,  conpute  M  from  Equation  1  and  then  L  from 
Equation  2  .  If  a  morning  phenomenon  (srnrise  or  the  beginning  of  morning 
twilight)  is  being  oonputed,  construct  an  initial  value  of  t  from  the  formula 

t  *  N  +  (6h  +  A  )/24 

vfriere  N  is  the  day  of  the  year*  and  A  is  the  observer's  longitude  expressed  in 
hours.  If  an  evening  phenomenon  is  being  computed,  use 

t  =  N  +  18h  +  A  )/24 

(b)  Step  2.  Solve  Equation  (3)  for  RA,  noting  that  RA  is  in  the  same  quadrant  at  L 
Transform  RA  to  hours  for  later  use  in  Equation  (6) 


*N0TE:  N  is  the  mmber  of  days  elapsed  since  0  January  of  the  current  year;  for  exanple,  N  *=  1  on  1 
January,  N  -  32  on  1  February,  and  so  forth.  Leap  year  days  are  counted,  so  N  -  366  on  31  Decenber  of  leap 
years  and  N  =  363  on  31  December  other  years.  N  is  listed  on  most  calendars.  If  desired,  it  can  be  com¬ 
puted  from: 


tot  (2Z5*5i) 


tot  (i4J&)  (1  +  tot 


')  +  Day  -  30 


where  Mon  is  month  (1  <:  Mxi  *  12),  Yr  is  the  year  (e.g. ,  1982),  and  Day  is  day  of 
month  (1  ^  Day  <  31). 
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(c)  Step  3.  Solve  Equation  4  for  sin  6  which  appears  in  Equation  5;  cos  6,  which 
also  is  required  in  Equation  5  should  be  determined  from  sin  6.  While  sin  5 
may  be  positive  or  negative,  cos  6  is  always  positive. 

(d)  Step  4.  Solve  Equation  5  for  H.  Since  computers  and  calculators  normally  give 
the  arccosine  in  the  range  0°  -  180°,  the  correct  quadrant  for  H  can  be  selected 
according  to: 

rising  phenomena,  H  =  360°  -  arccos  x; 
setting  phenomena,  H  =  arccos  x. 

In  other  words,  for  rising  phenomena,  H  trust  be  either  in  quadrant  3  or  5  (de¬ 
pending  on  the  sign  of  cos  H) ,  whereas  H  most  be  either  in  quadrant  1  or  2  for 
setting  phenomena.  Convert  H  from  degrees  to  hours  for  use  in  Equation  (6). 

(e)  Step  5.  Compute  T  from  Equation  6  ,  recalling  that  H  and  RA  must  be  expressed 
hours.  If  T  is  negative  or  greater  than  24*1,  it  should  be  converted  to  the  range 

_  24^  by  adding  or  subtracting  multiples  of  24^. 

(f)  Step  6.  Compute  UT  from  Equation  7  ,  where  X  must  be  expressed  in  hours.  UT  is 
an  approximation  to  the  time  of  sinrise,  sunset,  or  twilight,  referred  to  the 
Greenwich  meridian.  If  UT  is  greater  than  24*\  the  phenomenon  occurs  on  the 
following  day,  Greenwich  time.  If  UT  is  negative,  the  phenomenon  occurs  cn  the 
previous  day,  Greenwich  time. 

To  ensure  that  precision  is  not  lost  during  the  computations,  t  should  be  carried  to  four  decimal  places. 
Angles  should  be  expressed  to  three  decimals  of  a  degree  aid,  upon  conversion,  to  three  decimals  of  an 
hour.  Five  significant  digits  should  be  carried  for  the  trigonometric  functions .  Under  certain  conditions 
Equation  5  will  yield  a  value  of  |  cos  H  |  >  1,  indicating  the  absence  of  the  phenomenon  on  that  day. 

At  far  northern  latitudes,  for  example,  there  is  continuous  illunination  during  certain  summer  days  and 
continuous  darkness  during  winter  days.  EXAMPLE:  Compute  the  time  of  sunrise  on  25  June  at  Wayne,  New 
Jersey: 


Latitude:  40°. 9  North  Longitude:  74°. 3  West 

0  =  +49°. 9  sin  0  *  -H). 654744  cos  0  -  -K). 75585 

X  =  74. 3°/ 15  =  4.95h 

For  sunrise:  z  =  90°50'  cos  z  *=  -0.01454 

t  =  176d  +  (6h  +  4.95h)/24  =  176.456h 
M  -  0. 9356°(176.456d)  -  3.239°  =  170.626° 

L  =  170.626°  +  1.916° (0.16288)  +  0.020°  (-0.3214)  +  282.634° 

-  453.566°  -  93.566° 

d.  Solar  Coordinates.  The  true  geocentric  longitude  of  the  Srn  (L)  can  be  computed  to  an 
accuracy  of  +  1  minute  of  arc  from  the  following  formulas: 
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M  =  358°. 476  +  35999°. 050T 

L  =  279°. 691  +  36000°. 769T  +  (1°.919  -  0°.0048T)  sin  M  +  0°.020  sin  2M 

where  T  =  (N  -  2415010. 0)/36525  and  N  is  from  Equation  8-3. 

If  we  consider  the  Sun's  latitude  to  be  identically  zero,  the  right  ascension  (RA)  and  declination  (6)  of 
the  Sun  can  also  be  computed  to  +  1  minutes  of  arc  from 

tan  RA  =  cos  e  tan  L 

sin  6  =  sin  e  sin  L 

where  e,  the  obliquity  of  the  ecliptic,  can  be  computed  from  e  =  23.452°  -  0.013°T.  The  right  ascension 
is  always  in  the  same  quadrant  as  the  true  longitude.  Because  the  obliquity  varies  slowly,  a  single 
value  can  be  used  for  an  extended  period  of  time.  During  the  last  quarter  of  the  twentieth  century,  e  = 
23.441°  is  sufficiently  accurate.  Similarly  the  coefficient  of  sin  M  in  the  equation  for  L  changes 
slowly;  for  the  last  half  of  the  twentieth  century  a  value  of  1.916°  can  be  safety  used.  Although  there 
is  no  rigorous  limit  on  the  time  span  for  vhich  these  formulas  are  valid,  their  accuracy  gradually  deteri¬ 
orates  for  values  of  T  greater  than  a  couple  of  centuries. 

e.  Solar  Transit.  The  equation  of  time  (EqT)  is  the  hour  angle  of  the  true  Sun  minus  the  hour 
angle  of  the  mean  sun.  Thus,  it  is  the  difference:  apparent  solar  (sundial)  time  minus  meal  solar 
(clock)  time.  For  the  current  year  EqT  can  be  computed  to  an  accuracy  of  +  0.8  minutes  from  the  formula: 

(1)  EqT  =  -7.65msin  (0.9S56°t)  +0.53m  cos  (0.9356°t) 

-9.34m  sin  (1.9712°t)  -  2.92®  cos  (1.9712°t) 

where  t  is  the  number  of  days  since  0  January,  0^  UT.  If  higher  accuracy  is  required,  the  following 
formulas  will  give  EcfT  to  an  accuracy  of  +  2  seconds  during  the  current  year-. 

(2)  9  =  9.094°  +  0.98561°  t  +  1.916°  sin(0.9856°  t  -  3.562°) 

+  0.020°  sin  (1.9712°  t  -  7.124°) 

(3)  EqT  =  36.377m  +  3.94244m  t  -  4.CP  arctan  [(tan  9)/0.91747] 

where  t  is  the  number  of  days,  and  fractions  thereof,  since  0  January,  UT.  In  Equation  3  the  arc¬ 
tangent  should  yield  a  result  in  degrees  that  is  in  the  same  quadrant  as  0.  Near  the  end  of  the  year  0 
becomes  greater  than  360°.  When  this  occurs,  the  arctangent  in  Equation  3  should  also  be  greater  than 
360°.  Equations  2  and  3  can  be  used  to  compute  the  time  at  which  the  Sun  transits  the  local  meridian 
First  use  Equations  2  and  3  to  compute  EqT  for  t  =  N  +  (12^  +  x)/24,  where  N  is  the  day  of  the  year 
and  x  is  the  west  longitude  expressed  in  hours.  Then  the  local  mean  time  (HIT)  of  transit  is  given  to  an 
accuracy  of  +  2  seconds  by  LWT  =  12*1  -  EqT.  The  universal  time  of  local  transit  is  then  obtained  from 
UT  -  LMT  +  x-  EXAMPLE:  Compute  the  time  of  solar  transit  at  longitude  51°32.8^East  on  31  October  1983: 


X  “  51 . 547°/15  -  3.4365h  -  - 


For  solar  transit:  t  -  304d  +  (12h  -  3h.4365)/24  =  304.3568d 

9  =  9.094°  +  0.98561°  (304.3568d)  +  1.916°  (-0.8956) 

+  0.020°  (-0.7968)  =  307.339° 

EqT  =  36.377®  +  3/94244m  (304.3568d)  -  4.C P  arctan  [-1.21084/0.91747] 

=  36.377®  +  1199.908®  -  4.cP(304.989°)  =  +16.33® 

LMT  =  12h00®  -  16.33®  =  llh43.67® 

UT  =  llh43.67®  -  3h26.19®  =  8h17®29s 

f.  Moonrise  and  Moonset.  Time  of  moanrise  and  moonset  can  be  computed  for  specified  locations 
using  the  following  algorithm.  Between  latitudes  60°  North  and  60°  South,  the  phenomena  can  be  computed 
to  an  accuracy  of  +5®.  Although  the  algorithm  can  be  used  at  higher  latitudes,  its  accuracy  deteriorates 
near  dates  on  which  the  Mocn  remains  above  or  be  lew  the  horizon  for  more  than  twenty-four  hours. 

Notation:  0  =  latitude  of  observer  (north  is  positive;  south  is  negative); 

X  =  longitude  of  observer  (west  is  positive;  east  is  negative); 

■  i-th  approximation  to  universal  time  of  phenomenon,  expressed  in  days  from  0 
January,  0h  UT; 

GHA^  =  Moon's  GHA  at  time  t^;  'w 

5^  =  Moon’s  declination  at  time  t^  (north  is  positive;  south  is  negative); 

r^  =  i-th  correction  to  tQ,  thus  t^  =  tQ  +  r^; 

f+  =  i-th  approximation  to  Moon's  IHA  at  time  of  rise  or  set; 

4H.  =  i-th  approximation  to  Moon's  daily  rate  of  change  in  GHA. 

Fbrnaalas : 

(1)  4H.  *  (O^  =  GHAQ)/ri  for  i  -  0,  let  A  Hq  =  347.81° 

(2)  x^  *  cos  =  (.00233  -  sin  0  sin  6i)/(cos  0  cos  6^) 

(3)  r1+1  -  (Hi+1  -  Ho)/A  Ht 

<4>  ci+l  -  Co  +  Ci+1 

•  _ 

(a)  Step  1.  Let  t  =  N  +  (12^  +  x)/24,  vhere  N  is  the  day  of  the  year  and  x  is  the 
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Procedure: 


observer's  longitude  expressed  in  hours.  Set  i  =  0  and  begin  the 
following  iterative  process. 


(b)  Step  2.  For  time  t^  compute  the  Moan's  GHA  and  declination  to  navigational 
precision  (j<)'.l)  using  the  power  series  method  discussed  earlier.  Label  these 
quantities  GHA^  and  <5^,  respectively,  where  i  specifies  the  iteration  lumber. 
For  i  =  0,  compute  Hq  =  GHAq  -  x‘ 

(c)  Step  3.  If  i  =  0,  let  A  Hq  =  347.81°.  Otherwise  compute  A  from  Equation  1. 
If  A  Hi  <  0,  add  360°/  ^  to  A  Hi- 

(d)  Step  4.  Solve  Equation  2  for  H^+^.  Since  computers  and  calculators  normally 
give  the  arccosine  in  the  range  0°  -  180°,  the  correct  quadrant  for  LL+^  can  be 
selected  according  to: 


moonrise  computations,  =  360° 


arccos  x .+ ^ ; 


moonset  canputations,  =  arccos  x.+, . 


In  other  words,  near  the  time  of  moonrise  H._^  must  be  either  in  quadrant  3  or  quadrant  1  or  2.  For 
latitudes  higher  than  60°  (i.e.,  |0|>  60°),  the  condition  cos  |>|1  can  occur,  thereby  indicating  the 

absence  of  the  phenomenon  on  that  day. 


(e)  Step  5.  Compute  r^+1  from  Equation  3. 


If  r 


i+11 


^  rP1  ,  If  | r^ |  <  0.5  ,  proceed  to  Step  6. 

>0.5  ,  the  phenomenon  being  computed  occurs  on  the  day  before  the  day 


desired  (if  ri+^  is  negative)  or  on  the  day  following  the  day  desired  (if  r^  is 
positive) .  Normally,  the  phenomenon  on  the  desired  day  can  be  obtained  by  adding 
to  ri+^  (if  r^  is  negative),  or  subtracting  from  r^+^  (if  r^  is  positive), 

3 60° /ALL .  If  successful,  this  technique  will  produce  a  new  value  of  r^+^  in  the 
required  range.  However,  two  conditions  may  prevent  the  reduction  to  | r  |  < 


0.5 


d. 


for  low  values  of  i,r^+^  may  be  a  fairly  crude  approximation  to  the  ultimate 
value,  rn; 

each  month  there  is  one  day  (near  last  quarter)  on  which  there  is  no  moan- 
rise,  and  another  day  (near  first  quarter)  on  which  there  is  no  moonset. 

If  0.5d,  it  is  probably  worth  attempting  another  iteration  to  see  if  | r^ |  <  0.5d. 

(f)  Step  6.  Compute  t^+1  from  Equation  4.  If  |ti+1  -  tj  <  0.01d,  ti+1  is 

accurate  to  +  5ra.  Otherwise  it  is  necessary  to  iterate  the  solution  by  setting 
i  ■  i  +  1  and  executing  Steps  2  through  6  again. 


t  EXAMPLE:  Compute  moonrise  on  26  April  1983  at  latitude  38°59'  North  and  longitude  76°30'  West: 


0  =  4-  38.983  sin  0  =  40.62909  cos  0  =  40.77733 

X  =  +  76. 500°/ 15  =  4-  5.100^ 

The  day  is  found  to  be  day  116;  therefore, 

tQ  =  116d  4-  (12h  4-  5.100h)/24  =  116.71250d 

i  =  0:  Evaluating  the  power  series  for  tQl  using  constants  from  Figure  8-2 

GHAq  =  262.248°  <$0  =  -7.054° 

H  =  262.248°  -  76.500°  =  185.748° 
o 

A  HQ  =  347.81° 

xL  -  cos  Hq  =  [0.00233  -  (0.62909)  (-0.12280)  ]/[(0. 77733)  (0.99243)] 

-  40.10316  arccos  x^  =  84.079° 

Since  moonrise  is  sought,  is  in  quadrant  3  or  4; 

H.  =  275.921° 

rx  =  (275.921°  -  185. 748°)/ 347.81°  =  40.25926d 
|r^|  <  0.5das  required. 

tL  =  116. 71250d  4-  0. 25926d  =  116.97176d  =  26  April  23h19m  UT 
i  =  1:  Evaluating  the  pcwar  series  for  tp 
GHA1  =  352.576°  \  =  8.432° 

AHj^  =  (352.576°  -  262.248°)0.25926d  =  348.407° 

X2  =  cos  1^  =  [0.00233  -  (0.62909)(-0.14664)]/[(0.77733)(0. 98919) 

arccos  ^  =  82.935° 


=  40.12300 


Days  97  thru  102  JD  2445431 . 5  to  2445437 . 5  Dates  Apr  7  thru  Apr  12 

A  =  3.0  B  =  -33.33333333  W  =  97 


Moon 

Moan 

Moon 

M>cn 

GHA 

Dec 

H  P 

S  D 

Term 

o 

o 

O 

O 

0 

1293.8731 

-12.1168 

0.9163 

0.2497 

1 

1047.5575 

13.1999 

0.0252 

0.0069 

2 

0.7655 

2.5889 

0.0122 

0.0033 

3 

-0.6021 

-0.7075 

-0.0104 

-0.0028 

4 

-0.1375 

-0.0694 

-0.0048 

-0.0013 

5 

0.0302 

-0.0081 

0.0060 

0.0016 

Sums 

2341.7867 

2.8870 

0.9445 

0.2574 

Days  103  thru  108 

JD  2445437.5  to  2445443.5 

Dates  Apr  13  thru  Apr  18 

A  =  3.0 

B  =  -35.33333333 

W  =  103 

Moon 

Moan 

Moon 

Mb  on 

GHA 

Dec 

H  P 

S  D 

Term 

O 

O 

O 

O 

0 

1226.5320 

17.5393 

0.9703 

0.2644 

1 

1041.6475 

12.0323 

0.0227 

0.0062 

2 

-3.5078 

-4.0443 

-0.0078 

-0. 0021 

3 

-0.0457 

-1.4435 

-0.0023 

-0.0006 

4 

0.4431 

0.0674 

0.0023 

0.0006 

5 

0.0786 

0.0879 

0.0 

0.0 

Sums 

2265.1477 

24.2391 

0.9852 

0.2685 

Days  109  thru  114 

JD  2445443.5  to  2445449.5 

Dates  Apr  19  thru  Apr  24 

A  =  3.0 

B  =  -37.33333333 

W  =  109 

Moon 

Moon 

Moon 

Moan 

GHA 

Dec 

H  P 

S  D 

Term 

O 

O 

O 

O 

0 

1142.9034 

17.9258 

0.9882 

0.2692 

1 

1040.0961 

-12.1580 

-0.0032 

-0.0009 

2 

2.9739 

-4.6920 

-0.0073 

-0.0020 

3 

0.0904 

1.4035 

0.0004 

0.0001 

4 

-0.4803 

0.1619 

0.0002 

0.0 

5 

0.0652 

-0.0893 

-0.0012 

-0.0003 

Sums 

2185.6487 

2.5519 

0.9971 

0.2661 

Days  115  thru  120 

JD  2445449.5  to  2445455.5 

Dates  Apr  25  thru  Apr  30 

A  =  3.0 

B  =  -39.33333333 

W  =  115 

Moan 

Moon 

Moon 

Moan 

GHA 

Dec 

H  P 

S  D 

Term 

O 

O 

O 

O 

0 

1070.7377 

-13.5239 

0.9511 

0.2592 

1 

1044.7045 

-13.8267 

-0.0301 

-0.0082 

2 

-0.9260 

3.3560 

-0.0036 

-0.0010 

3 

-0.3224 

0.9767 

-0.0010 

-0.0003 

4 

0.2498 

-0.1510 

0.0015 

0.0004 

5 

0.0291 

-0.0206 

0.0034 

0.0009 

Suns 

2114.4727 

-23.1895 

0.9213 

0.2510 

Figure  8-2.  Pcwer  Series  Approximation  of  Nautical  Almanac  Data  for  Year  1983. 


Since  moonrise  is  sought,  H2  is  in  quadrant  3  or  4: 


^2  =  277.065° 

r2  =  (277/065°  -  185 . 748°) /348 . 407°  =  -K).26210d 

t2  =  116.71250d  +  0 . 26210d  =  116.97460d  =  26  April  23h26m  UT 

|t2  -  tL|  =  0.0028d  <  0.01d 

The  extremely  rapid  convergence  illustrated  in  this  example  occurs  frequently  but  not  invariably. 
Although  the  first  approximation  (t, )  will  often  give  adequate  precision  for  most  purposes,  it  is 

"*■  J 

reconmended  that  the  solution  be  iterated  and  that  the  convergence  criterion  -  tj  <0.01°)  be 

tested. 


9.  TOXIC  CORRIDOR  CALCULATIONS. 


The  theory  and  practice  of  toxic  corridor  calculations  are  discussed  at  length  in  AWS/TR-80/003,  and 
a  tutorial  refresher  on  the  topic  is  given  in  AWS/FM-81/007.  A  complete  program  for  the  TI-59  hand  calcu¬ 
lator  is  given  in  AWS/TR-80/003.  One  useful  model  is  presented  next,  for  illustration.  It  is  a  Non- 
Gaussian  dispersion  equation  extracted  from  AWS/FM-81/007.  Data  from  field  diffusion  studies  have  been 
used  not  only  to  develop  aid  evaluate  the  Gaussian  equation  but  also  to  develop  and  evaluate  statistical  ' 
equations  vbich  contain  those  meteorological  parameters  that  "explain"  or  fit  the  field  data  most  closely. 
Such  empirical  approaches  result  in  statistical  regression  equations.  Data  from  three  field  studies 
named  Ocean  Breeze,  Dry  Gulch,  and  Prairie  Grass  were  used  by  the  Air  Force  Cambridge  Research  Labora¬ 
tories  (new  the  Air  Force  Geophysics  Laboratory)  to  develop  an  equation  to  determine  downwind  peak 
concentration  of  airborne  contaminants  from  a  continuous  point  source.  This  empirically  derived  equation 
was  developed  from  data  collected  during  extensive  diffusion  experiments  with  tracer  released  simulating 
ground- level  continuous  point  sources.  Using  independent  data,  the  normalized  peak  concentrations  ob¬ 
tained  from  this  equation  have  been  found  to  be  accurate  within  a  factor  of  two,  65  percent  of  the  time 
and  within  a  factor  of  four,  94  percent  of  the  time.  The  equation  is 


Cp/Q  =1.75 


10 


_4x  -1-95  (aT  +  10)4’92 


If  one  is  concerned  with  the  downwind  distance,  X,  at  which  a  predetermined  concentration,  C  ,  will  occur 
for  a  known  source  strength,  Q,  and  temperature  difference,  delta-T  (A  T),  the  above  equation  can  be 
easily  solved  for  X.  This  has  been  done  below.  In  the  process,  appropriate  changes  were  made  to  the 
coefficient  to  convert  from  metric  units  to  English  units,  and  a  factor  was  added  to  convert  Cp/Q  from 
units  of  seconds  per  cubic  meter  to  units  of  PIW  per  lb/min.  The  converted  equation,  which  was  used  to 
generate  the  Toxic  Corridor  Length  Tables  in  AWS/TR-80/003,  is 


X 


3.28 


(  29.75) 

cam 


0.153  -0.513 

<^> 


-2.53  ~ 

(  A  T  +  10) 


where  X  =  Dowrwind  distance  in  feet.  As  used  here,  this  distance  defines  a  toxic  corridor  length.  This 
is  the  downwind  dimension  of  an  area  within  which  the  forecast  concentration  of  a  toxic  chemical  equals  or 
exceeds  a  specified  exposure  limit. 


P  =  A  probability  factor  used  to  determine  the  probability  that  a  specified 
concentration  is  not  exceeded  outside  the  corridor.  Calculations  in 
AWS/TR-30/003  assume  a  90-percent  probability;  therefore,  P  is  equal  to 
1.63.  Probability  factors  corresponding  to  other  probabilities  can  be 
found  in  Table  35,  AWS/TR-80/003. 

GM-J  =  Gram  molecular  weight  of  the  toxic  chemical. 

Cp  =  Peak  concentration  in  parts  per  million  by  volume  (PIW)  at  a  height  of 
approximately  5  feet  above  the  ground  at  a  given  downwind  distance,  X. 

Q  =  Source  strength  in  lb/min. 

A  T  =  The  tenperature  in  °F  at  54  feet  minus  the  temperature  at  6  feet. 

U0TE-.  A  negative  A  T  means  a  decrease  of  temperature  with  height  and 
a  positive  A  T  means  an  increase  with  height. 

In  this  equation,  the  meteorological  parameter  that  must  be  measured  or  estimated  is  temperature 
difference,  delta-T,  between  54  and  6  feet.  Delta-T  is  a  continuous  variable  that  may  be  approximated  as 
a  discrete  function  in  the  absence  of  delta-T  measuring  equipment.  Procedures  are  given  in  AWS/TR-80/003 
to  estimate  delta-T  if  delta-T  measuring  equipment  is  unavailable.  The  delta-T  input  to  this  equation 
should  not  be  confused  with  turbulence  typing  indexes  such  as  the  Pasquill-Gifford-Tumer  (PGT)  scheme 
used  in  AWS/TR-214.  The  PGT  index  is  used  to  estimate  the  stability  class  so  that  the  standard  deviations, 
o's,  of  contaminant  concentrations  can  be  estimated.  This  estimate  used  with  Gaussian  dispersion  models 
is  different  from  directly  in  putting  a  value  of  delta-T  into  the  Ocean  Breeze/Dry  Gulch  equation. 
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